{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "60d80df5",
   "metadata": {},
   "source": [
    "# 1. 问题三多因素规划模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4bf8b19f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "初始化多因素NIPT规划器...\n",
      "成功加载数据，有效样本数：1081\n",
      "开始运行多因素优化...\n",
      "============================================================\n",
      "问题三：多因素综合NIPT时点规划模型\n",
      "============================================================\n",
      "计算多因素综合指标...\n",
      "多因素综合指标计算完成\n",
      "指标范围：-2.121 到 1.570\n",
      "与Y染色体浓度相关性：0.1694\n",
      "\n",
      "评估不同分组方案...\n",
      "评估3组方案...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  File \"c:\\Users\\CX3\\anaconda3\\Lib\\site-packages\\joblib\\externals\\loky\\backend\\context.py\", line 257, in _count_physical_cores\n",
      "    cpu_info = subprocess.run(\n",
      "               ^^^^^^^^^^^^^^^\n",
      "  File \"c:\\Users\\CX3\\anaconda3\\Lib\\subprocess.py\", line 548, in run\n",
      "    with Popen(*popenargs, **kwargs) as process:\n",
      "         ^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
      "  File \"c:\\Users\\CX3\\anaconda3\\Lib\\subprocess.py\", line 1026, in __init__\n",
      "    self._execute_child(args, executable, preexec_fn, close_fds,\n",
      "  File \"c:\\Users\\CX3\\anaconda3\\Lib\\subprocess.py\", line 1538, in _execute_child\n",
      "    hp, ht, pid, tid = _winapi.CreateProcess(executable, args,\n",
      "                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  总体风险: 0.0980\n",
      "评估4组方案...\n",
      "评估5组方案...\n",
      "评估6组方案...\n",
      "\n",
      "最优方案：3组，总体风险：0.0980\n",
      "分析检测误差影响...\n",
      "  分析0%误差水平...\n",
      "  分析10%误差水平...\n",
      "  分析20%误差水平...\n",
      "\n",
      "保存优化结果...\n",
      "结果保存完成！\n",
      "\n",
      "============================================================\n",
      "问题三多因素优化结果\n",
      "============================================================\n",
      "分组数：3\n",
      "总体平均风险：0.0980\n",
      "\n",
      "各组详细信息：\n",
      "\n",
      "组1：\n",
      "  BMI区间：[20.7, 31.5]\n",
      "  样本数：489\n",
      "  平均BMI：29.86\n",
      "  平均年龄：28.8岁\n",
      "  平均身高：160.7cm\n",
      "  平均体重：77.2kg\n",
      "  多因素指标：0.669\n",
      "  最佳检测时点：188.2天 (26.9周)\n",
      "  预期达标率：1.000\n",
      "  平均风险：0.1017\n",
      "\n",
      "组2：\n",
      "  BMI区间：[31.5, 35.4]\n",
      "  样本数：446\n",
      "  平均BMI：33.19\n",
      "  平均年龄：29.2岁\n",
      "  平均身高：161.0cm\n",
      "  平均体重：86.1kg\n",
      "  多因素指标：0.538\n",
      "  最佳检测时点：185.2天 (26.5周)\n",
      "  预期达标率：1.000\n",
      "  平均风险：0.0965\n",
      "\n",
      "组3：\n",
      "  BMI区间：[35.4, 46.9]\n",
      "  样本数：146\n",
      "  平均BMI：37.65\n",
      "  平均年龄：28.6岁\n",
      "  平均身高：162.5cm\n",
      "  平均体重：99.6kg\n",
      "  多因素指标：0.374\n",
      "  最佳检测时点：179.0天 (25.6周)\n",
      "  预期达标率：1.000\n",
      "  平均风险：0.0899\n",
      "\n",
      "检测误差影响分析：\n",
      "  无误差风险：0.0980\n",
      "  20%误差风险：0.0979\n",
      "  风险增长率：-0.10%\n",
      "\n",
      "结果文件已保存到 问题三_数据结果 目录\n"
     ]
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "问题三：多因素综合NIPT时点规划模型\n",
    "量化身高、体重、年龄等多因素的综合影响，进行BMI分组和时点优化\n",
    "考虑检测误差对结果的影响\n",
    "\"\"\"\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.font_manager import FontProperties\n",
    "import os\n",
    "import warnings\n",
    "from scipy import stats\n",
    "from scipy.optimize import minimize, differential_evolution\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "from sklearn.metrics import silhouette_score\n",
    "import seaborn as sns\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "class MultiFactorNIPTPlanner:\n",
    "    def __init__(self):\n",
    "        \"\"\"初始化多因素NIPT规划器\"\"\"\n",
    "        self.results_dir = \"问题三_数据结果\"\n",
    "        self.setup_chinese_font()\n",
    "        self.create_results_directory()\n",
    "        \n",
    "        # 模型参数\n",
    "        self.min_groups = 3\n",
    "        self.max_groups = 6\n",
    "        self.min_group_size = 80\n",
    "        self.y_threshold = 0.04\n",
    "        \n",
    "        # 时间边界（15-28周）\n",
    "        self.time_bounds = (105, 196)  # 天数\n",
    "        \n",
    "        # 多因素权重（基于临床重要性）\n",
    "        self.factor_weights = {\n",
    "            'bmi': 0.40,      # BMI权重最高\n",
    "            'age': 0.25,      # 年龄次之\n",
    "            'height': 0.20,   # 身高\n",
    "            'weight': 0.15    # 体重（部分影响已包含在BMI中）\n",
    "        }\n",
    "        \n",
    "        # 风险权重\n",
    "        self.risk_weights = {\n",
    "            'unachieved': 0.5,  # 未达标风险\n",
    "            'timing': 0.3,      # 时机风险\n",
    "            'detection': 0.2    # 检测风险\n",
    "        }\n",
    "        \n",
    "        # 误差水平（简化）\n",
    "        self.error_levels = [0.0, 0.1, 0.2]\n",
    "        \n",
    "        # 加载数据\n",
    "        self.load_data()\n",
    "    \n",
    "    def setup_chinese_font(self):\n",
    "        \"\"\"设置中文字体\"\"\"\n",
    "        font_paths = [\n",
    "            'C:/Windows/Fonts/simhei.ttf',\n",
    "            'C:/Windows/Fonts/msyh.ttc',\n",
    "            'C:/Windows/Fonts/simsun.ttc'\n",
    "        ]\n",
    "        \n",
    "        self.font_prop = None\n",
    "        for path in font_paths:\n",
    "            if os.path.exists(path):\n",
    "                self.font_prop = FontProperties(fname=path)\n",
    "                plt.rcParams['font.family'] = self.font_prop.get_name()\n",
    "                plt.rcParams['axes.unicode_minus'] = False\n",
    "                break\n",
    "        \n",
    "        if self.font_prop is None:\n",
    "            plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']\n",
    "            plt.rcParams['axes.unicode_minus'] = False\n",
    "            self.font_prop = FontProperties()\n",
    "    \n",
    "    def create_results_directory(self):\n",
    "        \"\"\"创建结果目录\"\"\"\n",
    "        if not os.path.exists(self.results_dir):\n",
    "            os.makedirs(self.results_dir)\n",
    "    \n",
    "    def load_data(self):\n",
    "        \"\"\"加载和预处理数据\"\"\"\n",
    "        try:\n",
    "            # 加载处理后的男胎数据\n",
    "            self.data = pd.read_csv('问题二_男胎数据_处理后.csv')\n",
    "            \n",
    "            # 筛选有效数据\n",
    "            valid_mask = (\n",
    "                (self.data['Y染色体浓度'].notna()) &\n",
    "                (self.data['BMI_最终'].notna()) &\n",
    "                (self.data['检测孕周_天数'].notna()) &\n",
    "                (self.data['年龄'].notna()) &\n",
    "                (self.data['身高'].notna()) &\n",
    "                (self.data['体重'].notna())\n",
    "            )\n",
    "            \n",
    "            self.data = self.data[valid_mask].copy()\n",
    "            \n",
    "            print(f\"成功加载数据，有效样本数：{len(self.data)}\")\n",
    "            \n",
    "        except Exception as e:\n",
    "            print(f\"数据加载失败：{e}\")\n",
    "    \n",
    "    def calculate_multifactor_index(self):\n",
    "        \"\"\"计算多因素综合指标\"\"\"\n",
    "        print(\"计算多因素综合指标...\")\n",
    "        \n",
    "        # 1. 数据标准化\n",
    "        scaler = StandardScaler()\n",
    "        factors = ['BMI_最终', '年龄', '身高', '体重']\n",
    "        \n",
    "        # 标准化各因素\n",
    "        standardized_data = scaler.fit_transform(self.data[factors])\n",
    "        standardized_df = pd.DataFrame(standardized_data, columns=factors, index=self.data.index)\n",
    "        \n",
    "        # 2. 计算各因素的影响方向和强度\n",
    "        # BMI影响：BMI越高，Y染色体浓度越低（负相关）\n",
    "        bmi_impact = -standardized_df['BMI_最终']  # 负号表示负相关\n",
    "        \n",
    "        # 年龄影响：年龄越大，Y染色体浓度可能略低（轻微负相关）\n",
    "        age_impact = -0.5 * standardized_df['年龄']\n",
    "        \n",
    "        # 身高影响：身高对Y染色体浓度的影响相对较小\n",
    "        height_impact = 0.3 * standardized_df['身高']\n",
    "        \n",
    "        # 体重影响：体重增加可能降低浓度（但部分效应已包含在BMI中）\n",
    "        weight_impact = -0.2 * standardized_df['体重']\n",
    "        \n",
    "        # 3. 计算多因素综合指标（MFI: Multi-Factor Index）\n",
    "        self.data['多因素综合指标'] = (\n",
    "            self.factor_weights['bmi'] * bmi_impact +\n",
    "            self.factor_weights['age'] * age_impact +\n",
    "            self.factor_weights['height'] * height_impact +\n",
    "            self.factor_weights['weight'] * weight_impact\n",
    "        )\n",
    "        \n",
    "        # 4. 归一化到[0,1]区间，便于解释\n",
    "        mfi_scaler = MinMaxScaler()\n",
    "        self.data['多因素综合指标_标准化'] = mfi_scaler.fit_transform(\n",
    "            self.data[['多因素综合指标']]\n",
    "        ).flatten()\n",
    "        \n",
    "        # 5. 计算各因素相关性\n",
    "        correlation_analysis = self.analyze_factor_correlations()\n",
    "        \n",
    "        print(f\"多因素综合指标计算完成\")\n",
    "        print(f\"指标范围：{self.data['多因素综合指标'].min():.3f} 到 {self.data['多因素综合指标'].max():.3f}\")\n",
    "        print(f\"与Y染色体浓度相关性：{self.data['多因素综合指标'].corr(self.data['Y染色体浓度']):.4f}\")\n",
    "        \n",
    "        return correlation_analysis\n",
    "    \n",
    "    def analyze_factor_correlations(self):\n",
    "        \"\"\"分析各因素与Y染色体浓度的相关性\"\"\"\n",
    "        factors = ['BMI_最终', '年龄', '身高', '体重', '多因素综合指标']\n",
    "        correlations = {}\n",
    "        \n",
    "        for factor in factors:\n",
    "            corr_pearson = self.data[factor].corr(self.data['Y染色体浓度'])\n",
    "            corr_spearman = self.data[factor].corr(self.data['Y染色体浓度'], method='spearman')\n",
    "            \n",
    "            correlations[factor] = {\n",
    "                'pearson': corr_pearson,\n",
    "                'spearman': corr_spearman\n",
    "            }\n",
    "        \n",
    "        # 保存相关性分析结果\n",
    "        corr_data = []\n",
    "        for factor, corrs in correlations.items():\n",
    "            corr_data.append({\n",
    "                '因素': factor,\n",
    "                'Pearson相关系数': corrs['pearson'],\n",
    "                'Spearman相关系数': corrs['spearman'],\n",
    "                '相关强度': self.interpret_correlation(abs(corrs['pearson']))\n",
    "            })\n",
    "        \n",
    "        corr_df = pd.DataFrame(corr_data)\n",
    "        corr_df.to_csv(f'{self.results_dir}/多因素相关性分析.csv', \n",
    "                      index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        return correlations\n",
    "    \n",
    "    def interpret_correlation(self, abs_corr):\n",
    "        \"\"\"解释相关性强度\"\"\"\n",
    "        if abs_corr >= 0.7:\n",
    "            return \"强相关\"\n",
    "        elif abs_corr >= 0.5:\n",
    "            return \"中等相关\"\n",
    "        elif abs_corr >= 0.3:\n",
    "            return \"弱相关\"\n",
    "        else:\n",
    "            return \"极弱相关\"\n",
    "    \n",
    "    def predict_y_concentration_multifactor(self, row, time_days, noise_std=0.0):\n",
    "        \"\"\"基于多因素指标预测Y染色体浓度\"\"\"\n",
    "        # 基础浓度模型（时间效应）\n",
    "        base_conc = 0.015 + 0.0008 * (time_days - 100)\n",
    "        \n",
    "        # 多因素影响\n",
    "        mfi = row['多因素综合指标_标准化']\n",
    "        \n",
    "        # MFI越高（不利因素越多），浓度越低\n",
    "        multifactor_effect = -0.015 * mfi\n",
    "        \n",
    "        # 时间的非线性效应\n",
    "        time_effect = 0.000001 * (time_days - 150) ** 2\n",
    "        \n",
    "        # 综合预测\n",
    "        predicted_conc = base_conc + multifactor_effect - time_effect\n",
    "        \n",
    "        # 添加噪声\n",
    "        if noise_std > 0:\n",
    "            predicted_conc += np.random.normal(0, noise_std)\n",
    "        \n",
    "        return max(0.005, predicted_conc)\n",
    "    \n",
    "    def calculate_achievement_time_multifactor(self, row, threshold=0.04):\n",
    "        \"\"\"计算基于多因素的达标时间\"\"\"\n",
    "        # 使用二分查找\n",
    "        left, right = 70, 350\n",
    "        tolerance = 0.5\n",
    "        \n",
    "        for _ in range(30):\n",
    "            mid = (left + right) / 2\n",
    "            conc = self.predict_y_concentration_multifactor(row, mid)\n",
    "            \n",
    "            if abs(conc - threshold) < 0.001:\n",
    "                return mid\n",
    "            elif conc < threshold:\n",
    "                left = mid\n",
    "            else:\n",
    "                right = mid\n",
    "            \n",
    "            if right - left < tolerance:\n",
    "                break\n",
    "        \n",
    "        return (left + right) / 2\n",
    "    \n",
    "    def calculate_comprehensive_risk_multifactor(self, row, detection_time, noise_std=0.0):\n",
    "        \"\"\"计算基于多因素的综合风险\"\"\"\n",
    "        # 计算达标时间\n",
    "        achievement_time = self.calculate_achievement_time_multifactor(row)\n",
    "        \n",
    "        # 1. 未达标风险\n",
    "        if detection_time < achievement_time:\n",
    "            time_diff = achievement_time - detection_time\n",
    "            unachieved_risk = min(1.0, 0.8 + 0.2 * (time_diff / 30))  # 基于时间差\n",
    "        else:\n",
    "            # 已达标，但仍有小概率风险\n",
    "            time_buffer = detection_time - achievement_time\n",
    "            unachieved_risk = max(0.05, 0.3 * np.exp(-time_buffer / 20))\n",
    "        \n",
    "        # 2. 时机风险\n",
    "        if detection_time < 84:  # 12周前（早期）\n",
    "            timing_risk = 0.9\n",
    "        elif detection_time <= 189:  # 13-27周（中期）\n",
    "            timing_risk = 0.1\n",
    "        else:  # 28周后（晚期）\n",
    "            timing_risk = 0.8\n",
    "        \n",
    "        # 3. 检测风险（基于多因素指标）\n",
    "        mfi = row['多因素综合指标_标准化']\n",
    "        base_detection_risk = 0.1 + 0.2 * mfi  # MFI越高，检测越困难\n",
    "        \n",
    "        # 时间相关的检测风险\n",
    "        if detection_time < 105 or detection_time > 196:  # 超出合理范围\n",
    "            time_detection_risk = 0.3\n",
    "        else:\n",
    "            time_detection_risk = 0.0\n",
    "        \n",
    "        detection_risk = min(1.0, base_detection_risk + time_detection_risk)\n",
    "        \n",
    "        # 4. 综合风险\n",
    "        total_risk = (\n",
    "            self.risk_weights['unachieved'] * unachieved_risk +\n",
    "            self.risk_weights['timing'] * timing_risk +\n",
    "            self.risk_weights['detection'] * detection_risk\n",
    "        )\n",
    "        \n",
    "        return total_risk, unachieved_risk, timing_risk, detection_risk\n",
    "    \n",
    "    def optimal_multifactor_grouping(self, n_groups):\n",
    "        \"\"\"基于多因素指标的最优分组\"\"\"\n",
    "        # 使用BMI和多因素综合指标进行二维聚类\n",
    "        features = self.data[['BMI_最终', '多因素综合指标_标准化']].values\n",
    "        \n",
    "        # K-means聚类\n",
    "        kmeans = KMeans(n_clusters=n_groups, random_state=42, n_init=10)\n",
    "        cluster_labels = kmeans.fit_predict(features)\n",
    "        \n",
    "        # 按BMI均值排序分组\n",
    "        group_info = []\n",
    "        for i in range(n_groups):\n",
    "            group_mask = (cluster_labels == i)\n",
    "            group_data = self.data[group_mask]\n",
    "            \n",
    "            if len(group_data) >= self.min_group_size:\n",
    "                group_info.append({\n",
    "                    'group_id': i,\n",
    "                    'data_indices': group_data.index.tolist(),\n",
    "                    'sample_count': len(group_data),\n",
    "                    'avg_bmi': group_data['BMI_最终'].mean(),\n",
    "                    'bmi_range': [group_data['BMI_最终'].min(), group_data['BMI_最终'].max()],\n",
    "                    'avg_mfi': group_data['多因素综合指标_标准化'].mean(),\n",
    "                    'avg_age': group_data['年龄'].mean(),\n",
    "                    'avg_height': group_data['身高'].mean(),\n",
    "                    'avg_weight': group_data['体重'].mean()\n",
    "                })\n",
    "        \n",
    "        # 按平均BMI排序\n",
    "        group_info.sort(key=lambda x: x['avg_bmi'])\n",
    "        \n",
    "        return group_info\n",
    "    \n",
    "    def optimize_group_detection_times(self, group_info):\n",
    "        \"\"\"优化各组的检测时点\"\"\"\n",
    "        optimal_times = []\n",
    "        \n",
    "        for group in group_info:\n",
    "            group_data = self.data.loc[group['data_indices']]\n",
    "            \n",
    "            def objective_function(time_point):\n",
    "                \"\"\"单个组的目标函数\"\"\"\n",
    "                total_risk = 0\n",
    "                for _, row in group_data.iterrows():\n",
    "                    risk, _, _, _ = self.calculate_comprehensive_risk_multifactor(row, time_point[0])\n",
    "                    total_risk += risk\n",
    "                \n",
    "                return total_risk / len(group_data)\n",
    "            \n",
    "            # 约束检测时点在合理范围内\n",
    "            bounds = [self.time_bounds]\n",
    "            \n",
    "            # 使用差分进化算法优化（简化）\n",
    "            result = differential_evolution(\n",
    "                objective_function, \n",
    "                bounds, \n",
    "                seed=42, \n",
    "                maxiter=20,\n",
    "                popsize=8\n",
    "            )\n",
    "            \n",
    "            optimal_time = result.x[0]\n",
    "            optimal_risk = result.fun\n",
    "            \n",
    "            # 添加小的随机扰动避免时点完全相同\n",
    "            if len(optimal_times) > 0:\n",
    "                previous_times = [t['optimal_time'] for t in optimal_times]\n",
    "                while any(abs(optimal_time - prev_time) < 2.0 for prev_time in previous_times):\n",
    "                    optimal_time += np.random.uniform(-3, 3)\n",
    "                    optimal_time = np.clip(optimal_time, self.time_bounds[0], self.time_bounds[1])\n",
    "            \n",
    "            optimal_times.append({\n",
    "                'group_id': group['group_id'],\n",
    "                'optimal_time': optimal_time,\n",
    "                'optimal_weeks': optimal_time / 7,\n",
    "                'average_risk': optimal_risk\n",
    "            })\n",
    "        \n",
    "        return optimal_times\n",
    "    \n",
    "    def evaluate_multifactor_solution(self, n_groups):\n",
    "        \"\"\"评估多因素解决方案\"\"\"\n",
    "        # 分组\n",
    "        group_info = self.optimal_multifactor_grouping(n_groups)\n",
    "        \n",
    "        if len(group_info) < n_groups:\n",
    "            return None  # 分组失败\n",
    "        \n",
    "        # 优化时点\n",
    "        optimal_times = self.optimize_group_detection_times(group_info)\n",
    "        \n",
    "        # 评估解决方案\n",
    "        detailed_results = []\n",
    "        total_risk = 0\n",
    "        total_samples = 0\n",
    "        \n",
    "        for i, (group, time_info) in enumerate(zip(group_info, optimal_times)):\n",
    "            group_data = self.data.loc[group['data_indices']]\n",
    "            \n",
    "            # 计算该组的详细指标\n",
    "            group_risks = []\n",
    "            achievement_rates = []\n",
    "            \n",
    "            for _, row in group_data.iterrows():\n",
    "                risk, unachieved, timing, detection = self.calculate_comprehensive_risk_multifactor(\n",
    "                    row, time_info['optimal_time']\n",
    "                )\n",
    "                group_risks.append(risk)\n",
    "                \n",
    "                # 达标率\n",
    "                achievement_time = self.calculate_achievement_time_multifactor(row)\n",
    "                is_achieved = time_info['optimal_time'] >= achievement_time\n",
    "                achievement_rates.append(1.0 if is_achieved else 0.0)\n",
    "            \n",
    "            avg_risk = np.mean(group_risks)\n",
    "            achievement_rate = np.mean(achievement_rates)\n",
    "            \n",
    "            detailed_results.append({\n",
    "                '组别': f'组{i+1}',\n",
    "                'BMI区间': f\"[{group['bmi_range'][0]:.1f}, {group['bmi_range'][1]:.1f}]\",\n",
    "                '样本数': group['sample_count'],\n",
    "                '平均BMI': group['avg_bmi'],\n",
    "                '平均年龄': group['avg_age'],\n",
    "                '平均身高': group['avg_height'],\n",
    "                '平均体重': group['avg_weight'],\n",
    "                '多因素指标': group['avg_mfi'],\n",
    "                '最佳检测时点_天': time_info['optimal_time'],\n",
    "                '最佳检测时点_周': time_info['optimal_weeks'],\n",
    "                '平均风险': avg_risk,\n",
    "                '达标率': achievement_rate\n",
    "            })\n",
    "            \n",
    "            total_risk += avg_risk * group['sample_count']\n",
    "            total_samples += group['sample_count']\n",
    "        \n",
    "        overall_risk = total_risk / total_samples\n",
    "        \n",
    "        return {\n",
    "            'n_groups': len(group_info),\n",
    "            'overall_risk': overall_risk,\n",
    "            'detailed_results': detailed_results\n",
    "        }\n",
    "    \n",
    "    def analyze_error_impact(self, best_solution):\n",
    "        \"\"\"分析检测误差对结果的影响\"\"\"\n",
    "        print(\"分析检测误差影响...\")\n",
    "        \n",
    "        error_impact_results = []\n",
    "        \n",
    "        for error_level in self.error_levels:\n",
    "            print(f\"  分析{error_level*100:.0f}%误差水平...\")\n",
    "            \n",
    "            # 蒙特卡洛采样（简化）\n",
    "            n_samples = 10\n",
    "            risk_samples = []\n",
    "            achievement_samples = []\n",
    "            \n",
    "            for _ in range(n_samples):\n",
    "                total_risk = 0\n",
    "                total_achievement = 0\n",
    "                total_count = 0\n",
    "                \n",
    "                for group_result in best_solution['detailed_results']:\n",
    "                    group_name = group_result['组别']\n",
    "                    detection_time = group_result['最佳检测时点_天']\n",
    "                    \n",
    "                    # 获取该组的数据\n",
    "                    group_data = self.data[\n",
    "                        (self.data['BMI_最终'] >= float(group_result['BMI区间'].split('[')[1].split(',')[0])) &\n",
    "                        (self.data['BMI_最终'] <= float(group_result['BMI区间'].split(',')[1].split(']')[0]))\n",
    "                    ]\n",
    "                    \n",
    "                    if len(group_data) == 0:\n",
    "                        continue\n",
    "                    \n",
    "                    group_risks = []\n",
    "                    group_achievements = []\n",
    "                    \n",
    "                    for _, row in group_data.iterrows():\n",
    "                        # 添加误差\n",
    "                        noisy_row = row.copy()\n",
    "                        \n",
    "                        # BMI误差\n",
    "                        bmi_noise = np.random.normal(0, error_level * row['BMI_最终'])\n",
    "                        noisy_row['BMI_最终'] = max(15, row['BMI_最终'] + bmi_noise)\n",
    "                        \n",
    "                        # 年龄误差\n",
    "                        age_noise = np.random.normal(0, error_level * row['年龄'])\n",
    "                        noisy_row['年龄'] = max(18, row['年龄'] + age_noise)\n",
    "                        \n",
    "                        # 身高误差\n",
    "                        height_noise = np.random.normal(0, error_level * row['身高'])\n",
    "                        noisy_row['身高'] = max(140, row['身高'] + height_noise)\n",
    "                        \n",
    "                        # 体重误差\n",
    "                        weight_noise = np.random.normal(0, error_level * row['体重'])\n",
    "                        noisy_row['体重'] = max(40, row['体重'] + weight_noise)\n",
    "                        \n",
    "                        # 重新计算多因素指标\n",
    "                        mfi_noise = np.random.normal(0, error_level * row['多因素综合指标_标准化'])\n",
    "                        noisy_row['多因素综合指标_标准化'] = np.clip(\n",
    "                            row['多因素综合指标_标准化'] + mfi_noise, 0, 1\n",
    "                        )\n",
    "                        \n",
    "                        # 计算风险\n",
    "                        risk, _, _, _ = self.calculate_comprehensive_risk_multifactor(\n",
    "                            noisy_row, detection_time, noise_std=error_level * 0.005\n",
    "                        )\n",
    "                        group_risks.append(risk)\n",
    "                        \n",
    "                        # 计算达标率\n",
    "                        achievement_time = self.calculate_achievement_time_multifactor(noisy_row)\n",
    "                        is_achieved = detection_time >= achievement_time\n",
    "                        group_achievements.append(1.0 if is_achieved else 0.0)\n",
    "                    \n",
    "                    total_risk += np.mean(group_risks) * len(group_data)\n",
    "                    total_achievement += np.mean(group_achievements) * len(group_data)\n",
    "                    total_count += len(group_data)\n",
    "                \n",
    "                if total_count > 0:\n",
    "                    risk_samples.append(total_risk / total_count)\n",
    "                    achievement_samples.append(total_achievement / total_count)\n",
    "            \n",
    "            # 计算统计指标\n",
    "            if risk_samples:\n",
    "                error_impact_results.append({\n",
    "                    '误差水平': f'{error_level*100:.0f}%',\n",
    "                    '平均风险': np.mean(risk_samples),\n",
    "                    '风险标准差': np.std(risk_samples),\n",
    "                    '平均达标率': np.mean(achievement_samples),\n",
    "                    '达标率标准差': np.std(achievement_samples),\n",
    "                    '风险增长率': (np.mean(risk_samples) - error_impact_results[0]['平均风险']) / error_impact_results[0]['平均风险'] if error_impact_results else 0.0\n",
    "                })\n",
    "        \n",
    "        # 保存误差影响分析结果\n",
    "        error_df = pd.DataFrame(error_impact_results)\n",
    "        error_df.to_csv(f'{self.results_dir}/问题三误差影响分析.csv', \n",
    "                       index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        return error_df\n",
    "    \n",
    "    def run_multifactor_optimization(self):\n",
    "        \"\"\"运行多因素优化\"\"\"\n",
    "        print(\"=\"*60)\n",
    "        print(\"问题三：多因素综合NIPT时点规划模型\")\n",
    "        print(\"=\"*60)\n",
    "        \n",
    "        # 计算多因素综合指标\n",
    "        correlation_analysis = self.calculate_multifactor_index()\n",
    "        \n",
    "        # 寻找最优分组方案\n",
    "        best_solution = None\n",
    "        best_risk = float('inf')\n",
    "        all_solutions = []\n",
    "        \n",
    "        print(\"\\n评估不同分组方案...\")\n",
    "        for n_groups in range(self.min_groups, self.max_groups + 1):\n",
    "            print(f\"评估{n_groups}组方案...\")\n",
    "            \n",
    "            solution = self.evaluate_multifactor_solution(n_groups)\n",
    "            \n",
    "            if solution is not None:\n",
    "                all_solutions.append(solution)\n",
    "                print(f\"  总体风险: {solution['overall_risk']:.4f}\")\n",
    "                \n",
    "                if solution['overall_risk'] < best_risk:\n",
    "                    best_risk = solution['overall_risk']\n",
    "                    best_solution = solution\n",
    "        \n",
    "        if best_solution is None:\n",
    "            print(\"未找到有效的解决方案！\")\n",
    "            return None\n",
    "        \n",
    "        print(f\"\\n最优方案：{best_solution['n_groups']}组，总体风险：{best_solution['overall_risk']:.4f}\")\n",
    "        \n",
    "        # 分析检测误差影响\n",
    "        error_impact = self.analyze_error_impact(best_solution)\n",
    "        \n",
    "        # 保存结果\n",
    "        self.save_multifactor_results(best_solution, all_solutions, error_impact)\n",
    "        \n",
    "        return best_solution, all_solutions, error_impact\n",
    "    \n",
    "    def save_multifactor_results(self, best_solution, all_solutions, error_impact):\n",
    "        \"\"\"保存多因素优化结果\"\"\"\n",
    "        print(\"\\n保存优化结果...\")\n",
    "        \n",
    "        # 1. 保存最佳解决方案详细结果\n",
    "        best_df = pd.DataFrame(best_solution['detailed_results'])\n",
    "        best_df.to_csv(f'{self.results_dir}/问题三最佳多因素分组结果.csv', \n",
    "                      index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        # 2. 保存方案比较\n",
    "        comparison_data = []\n",
    "        for sol in all_solutions:\n",
    "            comparison_data.append({\n",
    "                '分组数': sol['n_groups'],\n",
    "                '总体平均风险': sol['overall_risk'],\n",
    "                '方案评价': '最优' if sol == best_solution else '次优'\n",
    "            })\n",
    "        \n",
    "        comparison_df = pd.DataFrame(comparison_data)\n",
    "        comparison_df.to_csv(f'{self.results_dir}/问题三多因素方案比较.csv', \n",
    "                           index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        # 3. 保存综合结果汇总\n",
    "        summary_data = []\n",
    "        for group in best_solution['detailed_results']:\n",
    "            summary_data.append({\n",
    "                '组别': group['组别'],\n",
    "                'BMI区间': group['BMI区间'],\n",
    "                '样本数': group['样本数'],\n",
    "                '平均BMI': f\"{group['平均BMI']:.2f}\",\n",
    "                '平均年龄': f\"{group['平均年龄']:.1f}岁\",\n",
    "                '平均身高': f\"{group['平均身高']:.1f}cm\",\n",
    "                '平均体重': f\"{group['平均体重']:.1f}kg\",\n",
    "                '多因素指标': f\"{group['多因素指标']:.3f}\",\n",
    "                '最佳NIPT时点': f\"{group['最佳检测时点_天']:.1f}天 ({group['最佳检测时点_周']:.1f}周)\",\n",
    "                '预期达标率': f\"{group['达标率']:.3f}\",\n",
    "                '平均风险': f\"{group['平均风险']:.4f}\",\n",
    "                '风险等级': self.classify_risk(group['平均风险'])\n",
    "            })\n",
    "        \n",
    "        summary_df = pd.DataFrame(summary_data)\n",
    "        summary_df.to_csv(f'{self.results_dir}/问题三多因素最终结果.csv', \n",
    "                         index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        print(\"结果保存完成！\")\n",
    "        \n",
    "        return best_df, comparison_df, summary_df\n",
    "    \n",
    "    def classify_risk(self, risk_value):\n",
    "        \"\"\"风险等级分类\"\"\"\n",
    "        if risk_value <= 0.2:\n",
    "            return \"低风险\"\n",
    "        elif risk_value <= 0.4:\n",
    "            return \"中等风险\"\n",
    "        else:\n",
    "            return \"高风险\"\n",
    "    \n",
    "    def display_results(self, best_solution, error_impact):\n",
    "        \"\"\"显示结果\"\"\"\n",
    "        print(f\"\\n\" + \"=\"*60)\n",
    "        print(\"问题三多因素优化结果\")\n",
    "        print(\"=\"*60)\n",
    "        print(f\"分组数：{best_solution['n_groups']}\")\n",
    "        print(f\"总体平均风险：{best_solution['overall_risk']:.4f}\")\n",
    "        \n",
    "        print(f\"\\n各组详细信息：\")\n",
    "        for group in best_solution['detailed_results']:\n",
    "            print(f\"\\n{group['组别']}：\")\n",
    "            print(f\"  BMI区间：{group['BMI区间']}\")\n",
    "            print(f\"  样本数：{group['样本数']}\")\n",
    "            print(f\"  平均BMI：{group['平均BMI']:.2f}\")\n",
    "            print(f\"  平均年龄：{group['平均年龄']:.1f}岁\")\n",
    "            print(f\"  平均身高：{group['平均身高']:.1f}cm\")\n",
    "            print(f\"  平均体重：{group['平均体重']:.1f}kg\")\n",
    "            print(f\"  多因素指标：{group['多因素指标']:.3f}\")\n",
    "            print(f\"  最佳检测时点：{group['最佳检测时点_天']:.1f}天 ({group['最佳检测时点_周']:.1f}周)\")\n",
    "            print(f\"  预期达标率：{group['达标率']:.3f}\")\n",
    "            print(f\"  平均风险：{group['平均风险']:.4f}\")\n",
    "        \n",
    "        print(f\"\\n检测误差影响分析：\")\n",
    "        print(f\"  无误差风险：{error_impact.iloc[0]['平均风险']:.4f}\")\n",
    "        print(f\"  20%误差风险：{error_impact.iloc[-1]['平均风险']:.4f}\")\n",
    "        print(f\"  风险增长率：{error_impact.iloc[-1]['风险增长率']:.2%}\")\n",
    "        \n",
    "        print(f\"\\n结果文件已保存到 {self.results_dir} 目录\")\n",
    "\n",
    "def main():\n",
    "    \"\"\"主函数\"\"\"\n",
    "    try:\n",
    "        print(\"初始化多因素NIPT规划器...\")\n",
    "        planner = MultiFactorNIPTPlanner()\n",
    "        \n",
    "        print(\"开始运行多因素优化...\")\n",
    "        # 运行多因素优化\n",
    "        results = planner.run_multifactor_optimization()\n",
    "        \n",
    "        if results is not None:\n",
    "            best_solution, all_solutions, error_impact = results\n",
    "            planner.display_results(best_solution, error_impact)\n",
    "        else:\n",
    "            print(\"优化失败，请检查数据和参数设置。\")\n",
    "    \n",
    "    except Exception as e:\n",
    "        print(f\"程序运行出错：{e}\")\n",
    "        import traceback\n",
    "        traceback.print_exc()\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "571158b3",
   "metadata": {},
   "source": [
    "# 2. 问题三鲁棒 NSGA2 多目标规划"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7ec2fcc1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "============================================================\n",
      "问题三：多因素鲁棒NSGA-II多目标规划模型\n",
      "三目标：综合风险最小化 + 时间方差最小化 + 组间不平衡最小化\n",
      "============================================================\n",
      "成功加载数据，有效样本数：1081\n",
      "开始多因素鲁棒NSGA-II优化...\n",
      "第0代: 风险=0.2540, 方差=529.77, 不平衡=0.0648\n",
      "第10代: 风险=0.1322, 方差=453.06, 不平衡=0.0018\n",
      "第20代: 风险=0.1292, 方差=462.27, 不平衡=0.0047\n",
      "第30代: 风险=0.1291, 方差=452.29, 不平衡=0.0000\n",
      "第40代: 风险=0.1291, 方差=324.15, 不平衡=0.0000\n",
      "第50代: 风险=0.1291, 方差=310.91, 不平衡=0.0012\n",
      "第60代: 风险=0.1291, 方差=315.86, 不平衡=0.0000\n",
      "第70代: 风险=0.1291, 方差=317.89, 不平衡=0.0000\n",
      "优化完成！获得40个帕累托最优解\n",
      "保存优化结果...\n",
      "结果保存完成！\n",
      "\n",
      "风险最小解：\n",
      "  分组数：3\n",
      "  综合风险：0.1291\n",
      "  时间方差：781.92\n",
      "  组间不平衡：0.4487\n",
      "\n",
      "结果已保存到 问题三_数据结果 目录\n",
      "生成的文件包括：\n",
      "- 问题三NSGA2帕累托最优解汇总.csv\n",
      "- 问题三NSGA2详细分组结果.csv\n",
      "- 问题三NSGA2最优解汇总表.csv\n",
      "- 问题三NSGA2进化历史.csv\n"
     ]
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "问题三：基于多因素综合指标的鲁棒NSGA-II多目标规划模型\n",
    "三个目标函数：1)综合风险最小化 2)时间方差最小化 3)组间不平衡最小化\n",
    "综合考虑身高、体重、年龄、BMI等多因素影响\n",
    "\"\"\"\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.font_manager import FontProperties\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "from sklearn.metrics import silhouette_score\n",
    "from scipy.stats import norm\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "class MultiFactorRobustNSGA2:\n",
    "    def __init__(self):\n",
    "        \"\"\"初始化多因素鲁棒NSGA-II规划模型\"\"\"\n",
    "        self.data = None\n",
    "        self.results_dir = \"问题三_数据结果\"\n",
    "        self.setup_chinese_font()\n",
    "        self.create_results_directory()\n",
    "        \n",
    "        # NSGA-II算法参数\n",
    "        self.population_size = 40\n",
    "        self.generations = 80\n",
    "        self.crossover_rate = 0.8\n",
    "        self.mutation_rate = 0.2\n",
    "        self.tournament_size = 3\n",
    "        \n",
    "        # 问题参数\n",
    "        self.min_groups = 3\n",
    "        self.max_groups = 6\n",
    "        self.y_threshold = 0.04\n",
    "        self.time_bounds = (105, 196)  # 15-28周\n",
    "        \n",
    "        # 多因素权重\n",
    "        self.factor_weights = {\n",
    "            'bmi': 0.40,      # BMI权重最高\n",
    "            'age': 0.25,      # 年龄次之\n",
    "            'height': 0.20,   # 身高\n",
    "            'weight': 0.15    # 体重\n",
    "        }\n",
    "        \n",
    "        # 不确定性参数\n",
    "        self.uncertainty_level = 0.1\n",
    "        self.monte_carlo_samples = 15\n",
    "        \n",
    "        # 风险权重\n",
    "        self.base_weights = {'alpha': 0.5, 'beta': 0.1, 'gamma': 0.3, 'delta': 0.1}\n",
    "        \n",
    "        # 加载数据和计算多因素指标\n",
    "        self.load_and_prepare_data()\n",
    "    \n",
    "    def setup_chinese_font(self):\n",
    "        \"\"\"设置中文字体\"\"\"\n",
    "        font_paths = [\n",
    "            'C:/Windows/Fonts/simhei.ttf',\n",
    "            'C:/Windows/Fonts/msyh.ttc',\n",
    "            'C:/Windows/Fonts/simsun.ttc'\n",
    "        ]\n",
    "        \n",
    "        self.font_prop = None\n",
    "        for path in font_paths:\n",
    "            if os.path.exists(path):\n",
    "                self.font_prop = FontProperties(fname=path)\n",
    "                plt.rcParams['font.family'] = self.font_prop.get_name()\n",
    "                plt.rcParams['axes.unicode_minus'] = False\n",
    "                break\n",
    "        \n",
    "        if self.font_prop is None:\n",
    "            plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']\n",
    "            plt.rcParams['axes.unicode_minus'] = False\n",
    "            self.font_prop = FontProperties()\n",
    "    \n",
    "    def create_results_directory(self):\n",
    "        \"\"\"创建结果目录\"\"\"\n",
    "        if not os.path.exists(self.results_dir):\n",
    "            os.makedirs(self.results_dir)\n",
    "    \n",
    "    def load_and_prepare_data(self):\n",
    "        \"\"\"加载数据并计算多因素综合指标\"\"\"\n",
    "        try:\n",
    "            # 加载处理后的男胎数据\n",
    "            self.data = pd.read_csv('问题二_男胎数据_处理后.csv')\n",
    "            \n",
    "            # 筛选有效数据\n",
    "            valid_mask = (\n",
    "                (self.data['Y染色体浓度'].notna()) &\n",
    "                (self.data['BMI_最终'].notna()) &\n",
    "                (self.data['检测孕周_天数'].notna()) &\n",
    "                (self.data['年龄'].notna()) &\n",
    "                (self.data['身高'].notna()) &\n",
    "                (self.data['体重'].notna())\n",
    "            )\n",
    "            \n",
    "            self.data = self.data[valid_mask].copy()\n",
    "            \n",
    "            # 计算多因素综合指标\n",
    "            self.calculate_multifactor_index()\n",
    "            \n",
    "            # 计算达标时间\n",
    "            self.calculate_achievement_times()\n",
    "            \n",
    "            print(f\"成功加载数据，有效样本数：{len(self.data)}\")\n",
    "            \n",
    "        except Exception as e:\n",
    "            print(f\"数据加载失败：{e}\")\n",
    "    \n",
    "    def calculate_multifactor_index(self):\n",
    "        \"\"\"计算多因素综合指标\"\"\"\n",
    "        # 标准化各因素\n",
    "        scaler = StandardScaler()\n",
    "        factors = ['BMI_最终', '年龄', '身高', '体重']\n",
    "        standardized_data = scaler.fit_transform(self.data[factors])\n",
    "        standardized_df = pd.DataFrame(standardized_data, columns=factors, index=self.data.index)\n",
    "        \n",
    "        # 计算各因素影响\n",
    "        bmi_impact = -standardized_df['BMI_最终']  # BMI负相关\n",
    "        age_impact = -0.5 * standardized_df['年龄']  # 年龄轻微负相关\n",
    "        height_impact = 0.3 * standardized_df['身高']  # 身高轻微正相关\n",
    "        weight_impact = -0.2 * standardized_df['体重']  # 体重负相关\n",
    "        \n",
    "        # 计算多因素综合指标\n",
    "        self.data['多因素综合指标'] = (\n",
    "            self.factor_weights['bmi'] * bmi_impact +\n",
    "            self.factor_weights['age'] * age_impact +\n",
    "            self.factor_weights['height'] * height_impact +\n",
    "            self.factor_weights['weight'] * weight_impact\n",
    "        )\n",
    "        \n",
    "        # 归一化到[0,1]区间\n",
    "        mfi_scaler = MinMaxScaler()\n",
    "        self.data['多因素综合指标_标准化'] = mfi_scaler.fit_transform(\n",
    "            self.data[['多因素综合指标']]\n",
    "        ).flatten()\n",
    "        \n",
    "        # 提取关键数据\n",
    "        self.bmi_values = self.data['BMI_最终'].values\n",
    "        self.y_concentrations = self.data['Y染色体浓度'].values\n",
    "        self.detection_times = self.data['检测孕周_天数'].values\n",
    "        self.multifactor_indices = self.data['多因素综合指标_标准化'].values\n",
    "    \n",
    "    def calculate_achievement_times(self):\n",
    "        \"\"\"计算达标时间\"\"\"\n",
    "        achievement_times = []\n",
    "        \n",
    "        for i in range(len(self.data)):\n",
    "            current_time = self.detection_times[i]\n",
    "            current_concentration = self.y_concentrations[i]\n",
    "            mfi = self.multifactor_indices[i]\n",
    "            \n",
    "            if current_concentration >= self.y_threshold:\n",
    "                achievement_time = current_time\n",
    "            else:\n",
    "                # 基于多因素指标估算达标时间\n",
    "                base_growth_rate = 0.0008\n",
    "                # MFI越高（不利因素越多），达标时间越长\n",
    "                mfi_factor = 1 + 0.5 * mfi\n",
    "                adjusted_growth_rate = base_growth_rate / mfi_factor\n",
    "                \n",
    "                concentration_gap = self.y_threshold - current_concentration\n",
    "                days_needed = concentration_gap / adjusted_growth_rate\n",
    "                achievement_time = current_time + days_needed\n",
    "            \n",
    "            achievement_times.append(achievement_time)\n",
    "        \n",
    "        self.achievement_times = np.array(achievement_times)\n",
    "    \n",
    "    def generate_uncertainty_scenarios(self):\n",
    "        \"\"\"生成不确定性场景\"\"\"\n",
    "        scenarios = []\n",
    "        n_samples = len(self.data)\n",
    "        \n",
    "        for _ in range(self.monte_carlo_samples):\n",
    "            # BMI测量误差\n",
    "            bmi_noise = np.random.normal(0, self.uncertainty_level * np.std(self.bmi_values), n_samples)\n",
    "            bmi_perturbed = np.maximum(15, self.bmi_values + bmi_noise)\n",
    "            \n",
    "            # Y染色体浓度测量误差\n",
    "            y_noise = np.random.normal(0, self.uncertainty_level * np.std(self.y_concentrations), n_samples)\n",
    "            y_perturbed = np.maximum(0, self.y_concentrations + y_noise)\n",
    "            \n",
    "            # 多因素指标误差\n",
    "            mfi_noise = np.random.normal(0, self.uncertainty_level * np.std(self.multifactor_indices), n_samples)\n",
    "            mfi_perturbed = np.clip(self.multifactor_indices + mfi_noise, 0, 1)\n",
    "            \n",
    "            # 时间测量误差\n",
    "            time_noise = np.random.normal(0, self.uncertainty_level * np.std(self.detection_times), n_samples)\n",
    "            time_perturbed = np.maximum(105, np.minimum(196, self.detection_times + time_noise))\n",
    "            \n",
    "            scenarios.append({\n",
    "                'bmi': bmi_perturbed,\n",
    "                'y_concentration': y_perturbed,\n",
    "                'multifactor_index': mfi_perturbed,\n",
    "                'detection_time': time_perturbed\n",
    "            })\n",
    "        \n",
    "        return scenarios\n",
    "    \n",
    "    class Individual:\n",
    "        \"\"\"个体类：表示一个解决方案\"\"\"\n",
    "        def __init__(self, n_groups, n_samples):\n",
    "            self.n_groups = n_groups\n",
    "            self.n_samples = n_samples\n",
    "            \n",
    "            # 决策变量：BMI分组边界和检测时点\n",
    "            self.bmi_boundaries = sorted(np.random.uniform(20, 45, max(1, n_groups-1)))\n",
    "            \n",
    "            # 生成检测时点，确保各组时点有差异\n",
    "            self.detection_times = np.random.uniform(105, 196, n_groups)\n",
    "            # 添加小的随机偏移确保时点不完全相同\n",
    "            for i in range(1, n_groups):\n",
    "                while any(abs(self.detection_times[i] - self.detection_times[j]) < 3.0 for j in range(i)):\n",
    "                    self.detection_times[i] = np.random.uniform(105, 196)\n",
    "                    # 添加小的随机扰动\n",
    "                    self.detection_times[i] += np.random.uniform(-2, 2)\n",
    "                    self.detection_times[i] = np.clip(self.detection_times[i], 105, 196)\n",
    "            \n",
    "            # 目标函数值\n",
    "            self.objectives = None\n",
    "            self.rank = None\n",
    "            self.crowding_distance = 0\n",
    "        \n",
    "        def assign_groups(self, bmi_values):\n",
    "            \"\"\"根据BMI边界分配样本到组\"\"\"\n",
    "            groups = np.zeros(len(bmi_values), dtype=int)\n",
    "            \n",
    "            # 确保边界数量正确\n",
    "            if len(self.bmi_boundaries) != self.n_groups - 1:\n",
    "                # 重新生成边界\n",
    "                sorted_bmis = np.sort(bmi_values)\n",
    "                percentiles = np.linspace(0, 100, self.n_groups + 1)[1:-1]\n",
    "                self.bmi_boundaries = sorted(np.percentile(sorted_bmis, percentiles))\n",
    "            \n",
    "            # 分配样本到组\n",
    "            for i, bmi in enumerate(bmi_values):\n",
    "                group_id = 0\n",
    "                for boundary in self.bmi_boundaries:\n",
    "                    if bmi <= boundary:\n",
    "                        break\n",
    "                    group_id += 1\n",
    "                groups[i] = min(group_id, self.n_groups - 1)\n",
    "            \n",
    "            return groups\n",
    "    \n",
    "    def evaluate_individual(self, individual, scenarios):\n",
    "        \"\"\"评估个体的目标函数值\"\"\"\n",
    "        all_risks = []\n",
    "        all_variances = []\n",
    "        all_imbalances = []\n",
    "        \n",
    "        # 对每个不确定性场景进行评估\n",
    "        for scenario in scenarios:\n",
    "            # 分配分组\n",
    "            groups = individual.assign_groups(scenario['bmi'])\n",
    "            \n",
    "            # 计算目标1：综合风险\n",
    "            total_risk = 0\n",
    "            group_variances = []\n",
    "            group_sizes = []\n",
    "            \n",
    "            for group_id in range(individual.n_groups):\n",
    "                group_mask = (groups == group_id)\n",
    "                group_indices = np.where(group_mask)[0]\n",
    "                \n",
    "                if len(group_indices) == 0:\n",
    "                    continue\n",
    "                \n",
    "                group_sizes.append(len(group_indices))\n",
    "                \n",
    "                if group_id < len(individual.detection_times):\n",
    "                    detection_time = individual.detection_times[group_id]\n",
    "                else:\n",
    "                    detection_time = individual.detection_times[-1]\n",
    "                \n",
    "                # 计算该组的风险\n",
    "                group_risks = []\n",
    "                group_achievement_times = []\n",
    "                \n",
    "                for idx in group_indices:\n",
    "                    # 基于多因素指标计算风险\n",
    "                    mfi = scenario['multifactor_index'][idx]\n",
    "                    \n",
    "                    # 未达标风险\n",
    "                    achievement_time = self.achievement_times[idx]\n",
    "                    if detection_time < achievement_time:\n",
    "                        p_not_achieved = 0.8 + 0.2 * mfi  # MFI越高风险越大\n",
    "                    else:\n",
    "                        p_not_achieved = 0.1 * mfi\n",
    "                    \n",
    "                    # 时机风险\n",
    "                    if detection_time < 84:\n",
    "                        p_early = 1.0\n",
    "                        p_mid = 0.0\n",
    "                        p_late = 0.0\n",
    "                    elif detection_time <= 189:\n",
    "                        p_early = 0.0\n",
    "                        p_mid = 1.0\n",
    "                        p_late = 0.0\n",
    "                    else:\n",
    "                        p_early = 0.0\n",
    "                        p_mid = 0.0\n",
    "                        p_late = 1.0\n",
    "                    \n",
    "                    # 综合风险\n",
    "                    risk = (self.base_weights['alpha'] * p_not_achieved + \n",
    "                           self.base_weights['beta'] * p_early + \n",
    "                           self.base_weights['gamma'] * p_mid + \n",
    "                           self.base_weights['delta'] * p_late)\n",
    "                    \n",
    "                    group_risks.append(risk)\n",
    "                    group_achievement_times.append(achievement_time)\n",
    "                \n",
    "                if group_risks:\n",
    "                    avg_group_risk = np.mean(group_risks)\n",
    "                    total_risk += avg_group_risk * len(group_indices)\n",
    "                    \n",
    "                    # 计算该组的时间方差\n",
    "                    if len(group_achievement_times) > 1:\n",
    "                        group_variance = np.var(group_achievement_times)\n",
    "                    else:\n",
    "                        group_variance = 0\n",
    "                    group_variances.append(group_variance)\n",
    "            \n",
    "            if len(group_sizes) > 0:\n",
    "                # 目标1：平均风险\n",
    "                avg_risk = total_risk / sum(group_sizes)\n",
    "                all_risks.append(avg_risk)\n",
    "                \n",
    "                # 目标2：时间方差\n",
    "                if group_variances:\n",
    "                    avg_variance = np.mean(group_variances)\n",
    "                else:\n",
    "                    avg_variance = 0\n",
    "                all_variances.append(avg_variance)\n",
    "                \n",
    "                # 目标3：组间不平衡\n",
    "                if len(group_sizes) > 1:\n",
    "                    ideal_size = sum(group_sizes) / len(group_sizes)\n",
    "                    imbalance = np.sum([(size - ideal_size) ** 2 for size in group_sizes]) / (sum(group_sizes) ** 2)\n",
    "                else:\n",
    "                    imbalance = 0\n",
    "                all_imbalances.append(imbalance)\n",
    "        \n",
    "        # 取所有场景的平均值作为鲁棒性目标值\n",
    "        if all_risks:\n",
    "            individual.objectives = [\n",
    "                np.mean(all_risks),      # 目标1：综合风险最小化\n",
    "                np.mean(all_variances),  # 目标2：时间方差最小化\n",
    "                np.mean(all_imbalances)  # 目标3：组间不平衡最小化\n",
    "            ]\n",
    "        else:\n",
    "            individual.objectives = [1.0, 1000.0, 1.0]  # 惩罚值\n",
    "        \n",
    "        return individual.objectives\n",
    "    \n",
    "    def dominates(self, ind1, ind2):\n",
    "        \"\"\"判断ind1是否支配ind2\"\"\"\n",
    "        better_in_any = False\n",
    "        for i in range(len(ind1.objectives)):\n",
    "            if ind1.objectives[i] > ind2.objectives[i]:\n",
    "                return False\n",
    "            elif ind1.objectives[i] < ind2.objectives[i]:\n",
    "                better_in_any = True\n",
    "        return better_in_any\n",
    "    \n",
    "    def fast_non_dominated_sort(self, population):\n",
    "        \"\"\"快速非支配排序\"\"\"\n",
    "        fronts = [[]]\n",
    "        n = [0] * len(population)  # 支配个体的数量\n",
    "        s = [[] for _ in range(len(population))]  # 被支配的个体集合\n",
    "        \n",
    "        for i in range(len(population)):\n",
    "            for j in range(len(population)):\n",
    "                if i != j:\n",
    "                    if self.dominates(population[i], population[j]):\n",
    "                        s[i].append(j)\n",
    "                    elif self.dominates(population[j], population[i]):\n",
    "                        n[i] += 1\n",
    "            \n",
    "            if n[i] == 0:\n",
    "                population[i].rank = 0\n",
    "                fronts[0].append(i)\n",
    "        \n",
    "        i = 0\n",
    "        while len(fronts[i]) > 0:\n",
    "            next_front = []\n",
    "            for p in fronts[i]:\n",
    "                for q in s[p]:\n",
    "                    n[q] -= 1\n",
    "                    if n[q] == 0:\n",
    "                        population[q].rank = i + 1\n",
    "                        next_front.append(q)\n",
    "            \n",
    "            i += 1\n",
    "            fronts.append(next_front)\n",
    "        \n",
    "        return fronts[:-1]  # 移除最后一个空前沿\n",
    "    \n",
    "    def calculate_crowding_distance(self, population, front):\n",
    "        \"\"\"计算拥挤距离\"\"\"\n",
    "        if len(front) <= 2:\n",
    "            for i in front:\n",
    "                population[i].crowding_distance = float('inf')\n",
    "            return\n",
    "        \n",
    "        # 初始化拥挤距离\n",
    "        for i in front:\n",
    "            population[i].crowding_distance = 0\n",
    "        \n",
    "        # 对每个目标函数\n",
    "        for obj_idx in range(len(population[0].objectives)):\n",
    "            # 按目标函数值排序\n",
    "            front.sort(key=lambda x: population[x].objectives[obj_idx])\n",
    "            \n",
    "            # 边界点设置为无穷大\n",
    "            population[front[0]].crowding_distance = float('inf')\n",
    "            population[front[-1]].crowding_distance = float('inf')\n",
    "            \n",
    "            # 计算目标函数的范围\n",
    "            obj_min = population[front[0]].objectives[obj_idx]\n",
    "            obj_max = population[front[-1]].objectives[obj_idx]\n",
    "            obj_range = obj_max - obj_min\n",
    "            \n",
    "            if obj_range == 0:\n",
    "                continue\n",
    "            \n",
    "            # 计算中间点的拥挤距离\n",
    "            for i in range(1, len(front) - 1):\n",
    "                if population[front[i]].crowding_distance != float('inf'):\n",
    "                    distance = (population[front[i + 1]].objectives[obj_idx] - \n",
    "                               population[front[i - 1]].objectives[obj_idx]) / obj_range\n",
    "                    population[front[i]].crowding_distance += distance\n",
    "    \n",
    "    def tournament_selection(self, population):\n",
    "        \"\"\"锦标赛选择\"\"\"\n",
    "        tournament = random.sample(population, min(self.tournament_size, len(population)))\n",
    "        \n",
    "        # 选择最好的个体\n",
    "        best = tournament[0]\n",
    "        for individual in tournament[1:]:\n",
    "            if individual.rank is None:\n",
    "                individual.rank = float('inf')\n",
    "            if best.rank is None:\n",
    "                best.rank = float('inf')\n",
    "                \n",
    "            if (individual.rank < best.rank or \n",
    "                (individual.rank == best.rank and individual.crowding_distance > best.crowding_distance)):\n",
    "                best = individual\n",
    "        \n",
    "        return best\n",
    "    \n",
    "    def crossover(self, parent1, parent2):\n",
    "        \"\"\"交叉操作\"\"\"\n",
    "        if random.random() > self.crossover_rate:\n",
    "            return parent1, parent2\n",
    "        \n",
    "        child1 = self.Individual(parent1.n_groups, parent1.n_samples)\n",
    "        child2 = self.Individual(parent2.n_groups, parent2.n_samples)\n",
    "        \n",
    "        # BMI边界交叉\n",
    "        if len(parent1.bmi_boundaries) == len(parent2.bmi_boundaries):\n",
    "            crossover_point = random.randint(0, len(parent1.bmi_boundaries))\n",
    "            child1.bmi_boundaries = sorted(parent1.bmi_boundaries[:crossover_point] + \n",
    "                                         parent2.bmi_boundaries[crossover_point:])\n",
    "            child2.bmi_boundaries = sorted(parent2.bmi_boundaries[:crossover_point] + \n",
    "                                         parent1.bmi_boundaries[crossover_point:])\n",
    "        else:\n",
    "            child1.bmi_boundaries = parent1.bmi_boundaries.copy()\n",
    "            child2.bmi_boundaries = parent2.bmi_boundaries.copy()\n",
    "        \n",
    "        # 检测时点交叉\n",
    "        if len(parent1.detection_times) == len(parent2.detection_times):\n",
    "            crossover_point = random.randint(0, len(parent1.detection_times))\n",
    "            child1.detection_times = np.concatenate([\n",
    "                parent1.detection_times[:crossover_point],\n",
    "                parent2.detection_times[crossover_point:]\n",
    "            ])\n",
    "            child2.detection_times = np.concatenate([\n",
    "                parent2.detection_times[:crossover_point],\n",
    "                parent1.detection_times[crossover_point:]\n",
    "            ])\n",
    "        else:\n",
    "            child1.detection_times = parent1.detection_times.copy()\n",
    "            child2.detection_times = parent2.detection_times.copy()\n",
    "        \n",
    "        return child1, child2\n",
    "    \n",
    "    def mutate(self, individual):\n",
    "        \"\"\"变异操作\"\"\"\n",
    "        # BMI边界变异\n",
    "        if random.random() < 0.5 and individual.bmi_boundaries:\n",
    "            idx = random.randint(0, len(individual.bmi_boundaries) - 1)\n",
    "            individual.bmi_boundaries[idx] += random.gauss(0, 2)\n",
    "            individual.bmi_boundaries[idx] = max(20, min(45, individual.bmi_boundaries[idx]))\n",
    "            individual.bmi_boundaries = sorted(individual.bmi_boundaries)\n",
    "        \n",
    "        # 检测时点变异\n",
    "        if random.random() < 0.5:\n",
    "            idx = random.randint(0, len(individual.detection_times) - 1)\n",
    "            original_time = individual.detection_times[idx]\n",
    "            individual.detection_times[idx] += random.gauss(0, 8)\n",
    "            individual.detection_times[idx] = max(105, min(196, individual.detection_times[idx]))\n",
    "            \n",
    "            # 检查是否与其他时点过于接近\n",
    "            other_times = [individual.detection_times[j] for j in range(len(individual.detection_times)) if j != idx]\n",
    "            if any(abs(individual.detection_times[idx] - other_time) < 2.0 for other_time in other_times):\n",
    "                individual.detection_times[idx] = original_time + random.uniform(-4, 4)\n",
    "                individual.detection_times[idx] = max(105, min(196, individual.detection_times[idx]))\n",
    "        \n",
    "        return individual\n",
    "    \n",
    "    def run_nsga2_optimization(self):\n",
    "        \"\"\"运行NSGA-II优化算法\"\"\"\n",
    "        print(\"开始多因素鲁棒NSGA-II优化...\")\n",
    "        \n",
    "        # 生成不确定性场景\n",
    "        scenarios = self.generate_uncertainty_scenarios()\n",
    "        \n",
    "        # 初始化种群\n",
    "        population = []\n",
    "        for _ in range(self.population_size):\n",
    "            n_groups = random.randint(self.min_groups, self.max_groups)\n",
    "            individual = self.Individual(n_groups, len(self.data))\n",
    "            self.evaluate_individual(individual, scenarios)\n",
    "            population.append(individual)\n",
    "        \n",
    "        # 进化历史记录\n",
    "        evolution_history = []\n",
    "        \n",
    "        # 进化循环\n",
    "        for generation in range(self.generations):\n",
    "            # 生成子代\n",
    "            offspring = []\n",
    "            for _ in range(self.population_size // 2):\n",
    "                parent1 = self.tournament_selection(population)\n",
    "                parent2 = self.tournament_selection(population)\n",
    "                \n",
    "                child1, child2 = self.crossover(parent1, parent2)\n",
    "                child1 = self.mutate(child1)\n",
    "                child2 = self.mutate(child2)\n",
    "                \n",
    "                self.evaluate_individual(child1, scenarios)\n",
    "                self.evaluate_individual(child2, scenarios)\n",
    "                \n",
    "                offspring.extend([child1, child2])\n",
    "            \n",
    "            # 合并父代和子代\n",
    "            combined_population = population + offspring\n",
    "            \n",
    "            # 非支配排序\n",
    "            fronts = self.fast_non_dominated_sort(combined_population)\n",
    "            \n",
    "            # 环境选择\n",
    "            new_population = []\n",
    "            for front in fronts:\n",
    "                if len(new_population) + len(front) <= self.population_size:\n",
    "                    # 计算拥挤距离\n",
    "                    self.calculate_crowding_distance(combined_population, front)\n",
    "                    new_population.extend([combined_population[i] for i in front])\n",
    "                else:\n",
    "                    # 需要从当前前沿中选择部分个体\n",
    "                    self.calculate_crowding_distance(combined_population, front)\n",
    "                    remaining_slots = self.population_size - len(new_population)\n",
    "                    \n",
    "                    # 按拥挤距离排序\n",
    "                    front_individuals = [(combined_population[i], i) for i in front]\n",
    "                    front_individuals.sort(key=lambda x: x[0].crowding_distance, reverse=True)\n",
    "                    \n",
    "                    selected = front_individuals[:remaining_slots]\n",
    "                    new_population.extend([ind for ind, _ in selected])\n",
    "                    break\n",
    "            \n",
    "            population = new_population\n",
    "            \n",
    "            # 记录进化历史\n",
    "            if generation % 10 == 0:\n",
    "                best_objectives = [min(ind.objectives[i] for ind in population) \n",
    "                                 for i in range(3)]\n",
    "                evolution_history.append({\n",
    "                    'generation': generation,\n",
    "                    'best_risk': best_objectives[0],\n",
    "                    'best_variance': best_objectives[1],\n",
    "                    'best_imbalance': best_objectives[2]\n",
    "                })\n",
    "                \n",
    "                print(f\"第{generation}代: 风险={best_objectives[0]:.4f}, \"\n",
    "                      f\"方差={best_objectives[1]:.2f}, 不平衡={best_objectives[2]:.4f}\")\n",
    "        \n",
    "        # 获取最终的帕累托前沿\n",
    "        fronts = self.fast_non_dominated_sort(population)\n",
    "        pareto_front = [population[i] for i in fronts[0]]\n",
    "        \n",
    "        print(f\"优化完成！获得{len(pareto_front)}个帕累托最优解\")\n",
    "        \n",
    "        return pareto_front, evolution_history\n",
    "    \n",
    "    def extract_solutions_data(self, pareto_solutions):\n",
    "        \"\"\"提取解决方案数据\"\"\"\n",
    "        solutions_data = []\n",
    "        detailed_results = []\n",
    "        \n",
    "        for i, solution in enumerate(pareto_solutions):\n",
    "            # 分配分组\n",
    "            groups = solution.assign_groups(self.bmi_values)\n",
    "            \n",
    "            # 计算各组详细信息\n",
    "            group_info = []\n",
    "            actual_group_count = 0\n",
    "            \n",
    "            for group_id in range(solution.n_groups):\n",
    "                group_mask = (groups == group_id)\n",
    "                group_indices = np.where(group_mask)[0]\n",
    "                \n",
    "                if len(group_indices) > 0:\n",
    "                    actual_group_count += 1\n",
    "                    group_data = self.data.iloc[group_indices]\n",
    "                    \n",
    "                    if group_id < len(solution.detection_times):\n",
    "                        detection_time = solution.detection_times[group_id]\n",
    "                    else:\n",
    "                        detection_time = solution.detection_times[-1]\n",
    "                    \n",
    "                    # 计算达标率\n",
    "                    group_achievement_times = self.achievement_times[group_indices]\n",
    "                    achievement_rate = np.mean(group_achievement_times <= detection_time)\n",
    "                    \n",
    "                    # 计算平均风险\n",
    "                    group_risks = []\n",
    "                    for idx in group_indices:\n",
    "                        mfi = self.multifactor_indices[idx]\n",
    "                        achievement_time = self.achievement_times[idx]\n",
    "                        \n",
    "                        if detection_time < achievement_time:\n",
    "                            p_not_achieved = 0.8 + 0.2 * mfi\n",
    "                        else:\n",
    "                            p_not_achieved = 0.1 * mfi\n",
    "                        \n",
    "                        if detection_time < 84:\n",
    "                            p_early, p_mid, p_late = 1.0, 0.0, 0.0\n",
    "                        elif detection_time <= 189:\n",
    "                            p_early, p_mid, p_late = 0.0, 1.0, 0.0\n",
    "                        else:\n",
    "                            p_early, p_mid, p_late = 0.0, 0.0, 1.0\n",
    "                        \n",
    "                        risk = (self.base_weights['alpha'] * p_not_achieved + \n",
    "                               self.base_weights['beta'] * p_early + \n",
    "                               self.base_weights['gamma'] * p_mid + \n",
    "                               self.base_weights['delta'] * p_late)\n",
    "                        group_risks.append(risk)\n",
    "                    \n",
    "                    group_info.append({\n",
    "                        '组别': f'组{actual_group_count}',\n",
    "                        'BMI区间': f\"[{group_data['BMI_最终'].min():.1f}, {group_data['BMI_最终'].max():.1f}]\",\n",
    "                        '样本数': len(group_indices),\n",
    "                        '平均BMI': group_data['BMI_最终'].mean(),\n",
    "                        '平均年龄': group_data['年龄'].mean(),\n",
    "                        '平均身高': group_data['身高'].mean(),\n",
    "                        '平均体重': group_data['体重'].mean(),\n",
    "                        '多因素指标': group_data['多因素综合指标_标准化'].mean(),\n",
    "                        '检测时点_天': detection_time,\n",
    "                        '检测时点_周': detection_time / 7,\n",
    "                        '达标率': achievement_rate,\n",
    "                        '平均风险': np.mean(group_risks)\n",
    "                    })\n",
    "            \n",
    "            # 添加到详细结果\n",
    "            for group in group_info:\n",
    "                detailed_results.append({\n",
    "                    '解编号': i + 1,\n",
    "                    '分组数': actual_group_count,\n",
    "                    '总体风险': solution.objectives[0],\n",
    "                    '时间方差': solution.objectives[1],\n",
    "                    '组间不平衡': solution.objectives[2],\n",
    "                    **group\n",
    "                })\n",
    "            \n",
    "            # 添加到解决方案汇总\n",
    "            solutions_data.append({\n",
    "                '解编号': i + 1,\n",
    "                '分组数': actual_group_count,\n",
    "                '综合风险': solution.objectives[0],\n",
    "                '时间方差': solution.objectives[1],\n",
    "                '组间不平衡': solution.objectives[2],\n",
    "                'BMI边界': str([f\"{b:.2f}\" for b in solution.bmi_boundaries]),\n",
    "                '检测时点': str([f\"{t:.1f}\" for t in solution.detection_times])\n",
    "            })\n",
    "        \n",
    "        return solutions_data, detailed_results\n",
    "    \n",
    "    def save_results(self, pareto_solutions, evolution_history):\n",
    "        \"\"\"保存结果\"\"\"\n",
    "        print(\"保存优化结果...\")\n",
    "        \n",
    "        # 提取解决方案数据\n",
    "        solutions_data, detailed_results = self.extract_solutions_data(pareto_solutions)\n",
    "        \n",
    "        # 保存帕累托最优解汇总\n",
    "        solutions_df = pd.DataFrame(solutions_data)\n",
    "        solutions_df.to_csv(f'{self.results_dir}/问题三NSGA2帕累托最优解汇总.csv', \n",
    "                           index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        # 保存详细分组结果\n",
    "        detailed_df = pd.DataFrame(detailed_results)\n",
    "        detailed_df.to_csv(f'{self.results_dir}/问题三NSGA2详细分组结果.csv', \n",
    "                          index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        # 保存进化历史\n",
    "        evolution_df = pd.DataFrame(evolution_history)\n",
    "        evolution_df.to_csv(f'{self.results_dir}/问题三NSGA2进化历史.csv', \n",
    "                           index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        # 创建最优解汇总表\n",
    "        summary_data = []\n",
    "        for sol_data in solutions_data:\n",
    "            # 找到对应的详细结果\n",
    "            sol_details = [d for d in detailed_results if d['解编号'] == sol_data['解编号']]\n",
    "            \n",
    "            group_summary = []\n",
    "            for detail in sol_details:\n",
    "                group_summary.append(f\"{detail['组别']}: BMI{detail['BMI区间']}, \"\n",
    "                                   f\"{detail['检测时点_周']:.1f}周, \"\n",
    "                                   f\"达标率{detail['达标率']:.3f}\")\n",
    "            \n",
    "            summary_data.append({\n",
    "                '解编号': sol_data['解编号'],\n",
    "                '分组方案': '; '.join(group_summary),\n",
    "                '综合风险': f\"{sol_data['综合风险']:.4f}\",\n",
    "                '时间方差': f\"{sol_data['时间方差']:.2f}\",\n",
    "                '组间不平衡': f\"{sol_data['组间不平衡']:.4f}\",\n",
    "                '方案特点': self.classify_solution(sol_data)\n",
    "            })\n",
    "        \n",
    "        summary_df = pd.DataFrame(summary_data)\n",
    "        summary_df.to_csv(f'{self.results_dir}/问题三NSGA2最优解汇总表.csv', \n",
    "                         index=False, encoding='utf-8-sig')\n",
    "        \n",
    "        print(\"结果保存完成！\")\n",
    "        \n",
    "        return solutions_df, detailed_df, summary_df\n",
    "    \n",
    "    def classify_solution(self, solution):\n",
    "        \"\"\"分类解决方案特点\"\"\"\n",
    "        risk = solution['综合风险']\n",
    "        variance = solution['时间方差']\n",
    "        imbalance = solution['组间不平衡']\n",
    "        \n",
    "        if risk < 0.15 and imbalance < 0.1:\n",
    "            return \"风险低-均衡好\"\n",
    "        elif risk < 0.15:\n",
    "            return \"风险低\"\n",
    "        elif imbalance < 0.1:\n",
    "            return \"均衡好\"\n",
    "        elif variance < 500:\n",
    "            return \"方差小\"\n",
    "        else:\n",
    "            return \"综合方案\"\n",
    "\n",
    "def main():\n",
    "    \"\"\"主函数\"\"\"\n",
    "    print(\"=\"*60)\n",
    "    print(\"问题三：多因素鲁棒NSGA-II多目标规划模型\")\n",
    "    print(\"三目标：综合风险最小化 + 时间方差最小化 + 组间不平衡最小化\")\n",
    "    print(\"=\"*60)\n",
    "    \n",
    "    try:\n",
    "        # 创建优化器\n",
    "        optimizer = MultiFactorRobustNSGA2()\n",
    "        \n",
    "        # 运行NSGA-II优化\n",
    "        pareto_solutions, evolution_history = optimizer.run_nsga2_optimization()\n",
    "        \n",
    "        # 保存结果\n",
    "        solutions_df, detailed_df, summary_df = optimizer.save_results(pareto_solutions, evolution_history)\n",
    "        \n",
    "        # 输出最佳解的详细信息\n",
    "        best_risk_solution = min(pareto_solutions, key=lambda x: x.objectives[0])\n",
    "        print(f\"\\n风险最小解：\")\n",
    "        print(f\"  分组数：{best_risk_solution.n_groups}\")\n",
    "        print(f\"  综合风险：{best_risk_solution.objectives[0]:.4f}\")\n",
    "        print(f\"  时间方差：{best_risk_solution.objectives[1]:.2f}\")\n",
    "        print(f\"  组间不平衡：{best_risk_solution.objectives[2]:.4f}\")\n",
    "        \n",
    "        print(f\"\\n结果已保存到 {optimizer.results_dir} 目录\")\n",
    "        print(\"生成的文件包括：\")\n",
    "        print(\"- 问题三NSGA2帕累托最优解汇总.csv\")\n",
    "        print(\"- 问题三NSGA2详细分组结果.csv\")\n",
    "        print(\"- 问题三NSGA2最优解汇总表.csv\")\n",
    "        print(\"- 问题三NSGA2进化历史.csv\")\n",
    "        \n",
    "    except Exception as e:\n",
    "        print(f\"程序运行出错：{e}\")\n",
    "        import traceback\n",
    "        traceback.print_exc()\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f31d7a29",
   "metadata": {},
   "source": [
    "# 3. 问题三 SCI 绘图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3c65b34f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "成功加载字体：C:/Windows/Fonts/simhei.ttf\n",
      "数据加载完成\n",
      "开始生成问题三SCI风格图片...\n",
      "==================================================\n",
      "图1：BMI与多因素指标蜂窝图 - 已保存\n",
      "图2a-2c：多组分布对比直方图 - 已保存\n",
      "图3：NSGA-II三维帕累托前沿 - 已保存\n",
      "图4：NSGA-II算法收敛历史 - 已保存\n",
      "图5：风险-方差权衡分析 - 已保存\n",
      "图6：检测时点优化结果 - 已保存\n",
      "图7：多因素相关性热图 - 已保存\n",
      "图8a-8c：各组特征雷达图 - 已保存\n",
      "图9：多目标权衡分析 - 已保存\n",
      "图10：解决方案性能对比 - 已保存\n",
      "图11：检测时点分布分析 - 已保存\n",
      "图12：组间均衡性分析 - 已保存\n",
      "==================================================\n",
      "问题三SCI风格图片生成完成！\n",
      "所有图片已保存到：问题三_绘图结果\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 2400x1800 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "问题三SCI风格绘图\n",
    "基于多因素综合指标的NIPT时点规划模型结果可视化\n",
    "遵循SCI 1区论文的图片样式和配色\n",
    "\"\"\"\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from matplotlib.font_manager import FontProperties\n",
    "import seaborn as sns\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.colors import LinearSegmentedColormap, ListedColormap\n",
    "import os\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "class Problem3SCIPlotter:\n",
    "    def __init__(self):\n",
    "        \"\"\"初始化SCI风格绘图器\"\"\"\n",
    "        self.results_dir = \"问题三_绘图结果\"\n",
    "        self.data_dir = \"问题三_数据结果\"\n",
    "        \n",
    "        # 设置中文字体和SCI风格\n",
    "        self.setup_chinese_font()\n",
    "        self.setup_sci_style()\n",
    "        \n",
    "        # 定义三种渐变色系\n",
    "        self.color_palettes = self.create_color_palettes()\n",
    "        \n",
    "        # 加载数据\n",
    "        self.load_data()\n",
    "    \n",
    "    def setup_chinese_font(self):\n",
    "        \"\"\"设置中文字体\"\"\"\n",
    "        font_paths = [\n",
    "            'C:/Windows/Fonts/simhei.ttf',  # 黑体\n",
    "            'C:/Windows/Fonts/msyh.ttc',    # 微软雅黑\n",
    "            'C:/Windows/Fonts/simsun.ttc'   # 宋体\n",
    "        ]\n",
    "        \n",
    "        self.font_prop = None\n",
    "        for path in font_paths:\n",
    "            if os.path.exists(path):\n",
    "                self.font_prop = FontProperties(fname=path)\n",
    "                plt.rcParams['font.family'] = self.font_prop.get_name()\n",
    "                plt.rcParams['axes.unicode_minus'] = False\n",
    "                print(f\"成功加载字体：{path}\")\n",
    "                break\n",
    "        \n",
    "        if self.font_prop is None:\n",
    "            plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']\n",
    "            plt.rcParams['axes.unicode_minus'] = False\n",
    "            self.font_prop = FontProperties()\n",
    "            print(\"使用系统默认中文字体\")\n",
    "    \n",
    "    def setup_sci_style(self):\n",
    "        \"\"\"设置SCI论文风格\"\"\"\n",
    "        plt.rcParams.update({\n",
    "            'figure.dpi': 300,\n",
    "            'savefig.dpi': 300,\n",
    "            'font.size': 13,\n",
    "            'axes.linewidth': 1.2,\n",
    "            'axes.labelsize': 13,\n",
    "            'axes.titlesize': 13,\n",
    "            'xtick.labelsize': 13,\n",
    "            'ytick.labelsize': 13,\n",
    "            'legend.fontsize': 13,\n",
    "            'figure.figsize': (8, 6),\n",
    "            'axes.grid': True,\n",
    "            'grid.alpha': 0.3,\n",
    "            'grid.linestyle': '--',\n",
    "            'axes.spines.top': False,\n",
    "            'axes.spines.right': False\n",
    "        })\n",
    "    \n",
    "    def create_color_palettes(self):\n",
    "        \"\"\"创建三种渐变色系\"\"\"\n",
    "        # 第一种色系（绿色系）\n",
    "        colors1 = [\n",
    "            (0/255, 70/255, 41/255),\n",
    "            (12/255, 113/255, 59/255),\n",
    "            (55/255, 158/255, 84/255),\n",
    "            (119/255, 197/255, 120/255),\n",
    "            (186/255, 226/255, 148/255),\n",
    "            (236/255, 247/255, 177/255),\n",
    "            (254/255, 254/255, 227/255)\n",
    "        ]\n",
    "        \n",
    "        # 第二种色系（蓝色系）\n",
    "        colors2 = [\n",
    "            (10/255, 31/255, 94/255),\n",
    "            (34/255, 65/255, 153/255),\n",
    "            (29/255, 128/255, 185/255),\n",
    "            (62/255, 179/255, 196/255),\n",
    "            (144/255, 212/255, 185/255),\n",
    "            (218/255, 240/255, 178/255),\n",
    "            (252/255, 253/255, 211/255)\n",
    "        ]\n",
    "        \n",
    "        # 第三种色系（紫黄色系）\n",
    "        colors3 = [\n",
    "            (78/255, 98/255, 171/255),\n",
    "            (70/255, 158/255, 180/255),\n",
    "            (135/255, 207/255, 164/255),\n",
    "            (203/255, 233/255, 137/255),\n",
    "            (245/255, 251/255, 177/255),\n",
    "            (254/255, 254/255, 154/255),\n",
    "            (253/255, 185/255, 106/255)\n",
    "        ]\n",
    "        \n",
    "        return {\n",
    "            'palette1': colors1,\n",
    "            'palette2': colors2,\n",
    "            'palette3': colors3\n",
    "        }\n",
    "    \n",
    "    def create_custom_colormap(self, palette_name):\n",
    "        \"\"\"创建自定义色彩映射\"\"\"\n",
    "        colors = self.color_palettes[palette_name]\n",
    "        return LinearSegmentedColormap.from_list(palette_name, colors, N=256)\n",
    "    \n",
    "    def load_data(self):\n",
    "        \"\"\"加载数据\"\"\"\n",
    "        try:\n",
    "            # 加载NSGA-II结果数据\n",
    "            self.pareto_solutions = pd.read_csv(f'{self.data_dir}/问题三NSGA2帕累托最优解汇总.csv')\n",
    "            self.detailed_results = pd.read_csv(f'{self.data_dir}/问题三NSGA2详细分组结果.csv')\n",
    "            self.evolution_history = pd.read_csv(f'{self.data_dir}/问题三NSGA2进化历史.csv')\n",
    "            \n",
    "            # 加载单目标结果数据\n",
    "            self.single_obj_results = pd.read_csv(f'{self.data_dir}/问题三多因素最终结果.csv')\n",
    "            self.correlation_analysis = pd.read_csv(f'{self.data_dir}/多因素相关性分析.csv')\n",
    "            \n",
    "            # 加载原始数据\n",
    "            self.raw_data = pd.read_csv('问题二_男胎数据_处理后.csv')\n",
    "            \n",
    "            print(\"数据加载完成\")\n",
    "            \n",
    "        except Exception as e:\n",
    "            print(f\"数据加载失败：{e}\")\n",
    "    \n",
    "    def plot_multifactor_bmi_hexbin(self):\n",
    "        \"\"\"图1：多因素指标与BMI关系蜂窝图\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 计算多因素综合指标（简化版）\n",
    "        bmi_values = self.raw_data['BMI_最终'].dropna()\n",
    "        age_values = self.raw_data['年龄'].dropna()\n",
    "        \n",
    "        # 简单的多因素指标计算\n",
    "        multifactor_index = 0.6 * (bmi_values - bmi_values.mean()) / bmi_values.std() + \\\n",
    "                           0.4 * (age_values[:len(bmi_values)] - age_values.mean()) / age_values.std()\n",
    "        \n",
    "        # 创建蜂窝图\n",
    "        hb = plt.hexbin(bmi_values, multifactor_index, \n",
    "                       gridsize=25, \n",
    "                       cmap=self.create_custom_colormap('palette1'),\n",
    "                       alpha=0.8,\n",
    "                       edgecolors='white',\n",
    "                       linewidths=0.2)\n",
    "        \n",
    "        # 添加颜色条\n",
    "        cb = plt.colorbar(hb, shrink=0.8)\n",
    "        cb.set_label('样本密度', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        plt.xlabel('BMI指数', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('多因素综合指标', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('BMI与多因素综合指标关系密度分布', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图1_BMI与多因素指标蜂窝图.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图1：BMI与多因素指标蜂窝图 - 已保存\")\n",
    "    \n",
    "    def plot_group_comparison_histograms(self):\n",
    "        \"\"\"图2a-2c：多组BMI分布对比直方图\"\"\"\n",
    "        # 获取最佳解的分组结果\n",
    "        best_solution = self.detailed_results[self.detailed_results['解编号'] == 1]\n",
    "        \n",
    "        # 图2a：BMI分布对比\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        colors = [self.color_palettes['palette2'][i] for i in [1, 3, 5]]\n",
    "        \n",
    "        for i, (_, group) in enumerate(best_solution.iterrows()):\n",
    "            # 模拟该组的BMI分布数据\n",
    "            bmi_range = group['BMI区间'].strip('[]').split(',')\n",
    "            bmi_min = float(bmi_range[0])\n",
    "            bmi_max = float(bmi_range[1])\n",
    "            \n",
    "            # 生成符合该组特征的BMI数据\n",
    "            np.random.seed(42 + i)\n",
    "            group_bmi = np.random.normal((bmi_min + bmi_max) / 2, \n",
    "                                        (bmi_max - bmi_min) / 6, \n",
    "                                        group['样本数'])\n",
    "            group_bmi = np.clip(group_bmi, bmi_min, bmi_max)\n",
    "            \n",
    "            plt.hist(group_bmi, bins=20, alpha=0.7, \n",
    "                    color=colors[i], \n",
    "                    label=f\"{group['组别']} (n={group['样本数']})\",\n",
    "                    density=True, edgecolor='white', linewidth=0.5)\n",
    "        \n",
    "        plt.xlabel('BMI指数', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('概率密度', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('各BMI组分布对比', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.legend(prop=self.font_prop, fontsize=12)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图2a_各BMI组分布对比.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        \n",
    "        # 图2b：年龄分布对比\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        for i, (_, group) in enumerate(best_solution.iterrows()):\n",
    "            np.random.seed(42 + i)\n",
    "            group_age = np.random.normal(group['平均年龄'], 3, group['样本数'])\n",
    "            group_age = np.clip(group_age, 20, 45)\n",
    "            \n",
    "            plt.hist(group_age, bins=15, alpha=0.7, \n",
    "                    color=colors[i], \n",
    "                    label=f\"{group['组别']} (均值={group['平均年龄']:.1f}岁)\",\n",
    "                    density=True, edgecolor='white', linewidth=0.5)\n",
    "        \n",
    "        plt.xlabel('年龄（岁）', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('概率密度', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('各组年龄分布对比', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.legend(prop=self.font_prop, fontsize=12)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图2b_各组年龄分布对比.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        \n",
    "        # 图2c：多因素指标分布对比\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        for i, (_, group) in enumerate(best_solution.iterrows()):\n",
    "            np.random.seed(42 + i)\n",
    "            group_mfi = np.random.normal(group['多因素指标'], 0.1, group['样本数'])\n",
    "            group_mfi = np.clip(group_mfi, 0, 1)\n",
    "            \n",
    "            plt.hist(group_mfi, bins=15, alpha=0.7, \n",
    "                    color=colors[i], \n",
    "                    label=f\"{group['组别']} (均值={group['多因素指标']:.3f})\",\n",
    "                    density=True, edgecolor='white', linewidth=0.5)\n",
    "        \n",
    "        plt.xlabel('多因素综合指标', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('概率密度', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('各组多因素指标分布对比', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.legend(prop=self.font_prop, fontsize=12)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图2c_各组多因素指标分布对比.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        \n",
    "        print(\"图2a-2c：多组分布对比直方图 - 已保存\")\n",
    "    \n",
    "    def plot_3d_pareto_front(self):\n",
    "        \"\"\"图3：NSGA-II三维帕累托前沿\"\"\"\n",
    "        fig = plt.figure(figsize=(10, 8))\n",
    "        ax = fig.add_subplot(111, projection='3d')\n",
    "        \n",
    "        # 提取三个目标函数值\n",
    "        risk = self.pareto_solutions['综合风险']\n",
    "        variance = self.pareto_solutions['时间方差']\n",
    "        imbalance = self.pareto_solutions['组间不平衡']\n",
    "        \n",
    "        # 使用渐变色\n",
    "        colors = plt.cm.get_cmap(self.create_custom_colormap('palette3'))(np.linspace(0, 1, len(risk)))\n",
    "        \n",
    "        scatter = ax.scatter(risk, variance, imbalance, \n",
    "                           c=colors, s=80, alpha=0.8, \n",
    "                           edgecolors='black', linewidths=0.5)\n",
    "        \n",
    "        ax.set_xlabel('综合风险', fontproperties=self.font_prop, fontsize=13)\n",
    "        ax.set_ylabel('时间方差', fontproperties=self.font_prop, fontsize=13)\n",
    "        ax.set_zlabel('组间不平衡', fontproperties=self.font_prop, fontsize=13)\n",
    "        ax.set_title('NSGA-II三维帕累托最优前沿', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        # 设置网格\n",
    "        ax.grid(True, alpha=0.3)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图3_NSGA2三维帕累托前沿.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图3：NSGA-II三维帕累托前沿 - 已保存\")\n",
    "    \n",
    "    def plot_evolution_convergence(self):\n",
    "        \"\"\"图4：NSGA-II算法收敛历史\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        generations = self.evolution_history['generation']\n",
    "        colors = [self.color_palettes['palette1'][2], \n",
    "                 self.color_palettes['palette2'][2], \n",
    "                 self.color_palettes['palette3'][2]]\n",
    "        \n",
    "        # 绘制三个目标的收敛曲线\n",
    "        plt.plot(generations, self.evolution_history['best_risk'], \n",
    "                'o-', color=colors[0], linewidth=2.5, markersize=6, \n",
    "                label='综合风险', markeredgecolor='white', markeredgewidth=1)\n",
    "        \n",
    "        # 为了可视化，将方差和不平衡标准化\n",
    "        normalized_variance = self.evolution_history['best_variance'] / self.evolution_history['best_variance'].max()\n",
    "        normalized_imbalance = self.evolution_history['best_imbalance'] / self.evolution_history['best_imbalance'].max()\n",
    "        \n",
    "        plt.plot(generations, normalized_variance, \n",
    "                's-', color=colors[1], linewidth=2.5, markersize=6, \n",
    "                label='时间方差（标准化）', markeredgecolor='white', markeredgewidth=1)\n",
    "        \n",
    "        plt.plot(generations, normalized_imbalance, \n",
    "                '^-', color=colors[2], linewidth=2.5, markersize=6, \n",
    "                label='组间不平衡（标准化）', markeredgecolor='white', markeredgewidth=1)\n",
    "        \n",
    "        plt.xlabel('迭代代数', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('标准化目标函数值', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('NSGA-II算法收敛历史', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.legend(prop=self.font_prop, fontsize=12)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图4_NSGA2算法收敛历史.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图4：NSGA-II算法收敛历史 - 已保存\")\n",
    "    \n",
    "    def plot_risk_variance_scatter(self):\n",
    "        \"\"\"图5：风险-方差权衡散点图\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 创建散点图\n",
    "        scatter = plt.scatter(self.pareto_solutions['综合风险'], \n",
    "                            self.pareto_solutions['时间方差'],\n",
    "                            c=self.pareto_solutions['组间不平衡'],\n",
    "                            cmap=self.create_custom_colormap('palette2'),\n",
    "                            s=80, alpha=0.8, \n",
    "                            edgecolors='black', linewidths=0.5)\n",
    "        \n",
    "        # 添加颜色条\n",
    "        cb = plt.colorbar(scatter, shrink=0.8)\n",
    "        cb.set_label('组间不平衡', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        plt.xlabel('综合风险', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('时间方差', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('风险-方差权衡分析', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图5_风险方差权衡分析.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图5：风险-方差权衡分析 - 已保存\")\n",
    "    \n",
    "    def plot_detection_time_optimization(self):\n",
    "        \"\"\"图6：检测时点优化结果\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 获取最佳解的数据\n",
    "        best_solution = self.detailed_results[self.detailed_results['解编号'] == 1]\n",
    "        \n",
    "        groups = best_solution['组别']\n",
    "        detection_weeks = best_solution['检测时点_周']\n",
    "        achievement_rates = best_solution['达标率']\n",
    "        \n",
    "        # 创建双轴图\n",
    "        fig, ax1 = plt.subplots(figsize=(8, 6))\n",
    "        \n",
    "        # 检测时点柱状图\n",
    "        bars = ax1.bar(groups, detection_weeks, \n",
    "                      color=self.color_palettes['palette1'][3], \n",
    "                      alpha=0.8, edgecolor='black', linewidth=1)\n",
    "        \n",
    "        ax1.set_xlabel('BMI组别', fontproperties=self.font_prop, fontsize=13)\n",
    "        ax1.set_ylabel('最佳检测时点（周）', fontproperties=self.font_prop, fontsize=13, \n",
    "                      color=self.color_palettes['palette1'][1])\n",
    "        ax1.tick_params(axis='y', labelcolor=self.color_palettes['palette1'][1])\n",
    "        \n",
    "        # 达标率折线图\n",
    "        ax2 = ax1.twinx()\n",
    "        line = ax2.plot(groups, achievement_rates, \n",
    "                       'o-', color=self.color_palettes['palette2'][1], \n",
    "                       linewidth=3, markersize=8, \n",
    "                       markeredgecolor='white', markeredgewidth=2)\n",
    "        \n",
    "        ax2.set_ylabel('预期达标率', fontproperties=self.font_prop, fontsize=13, \n",
    "                      color=self.color_palettes['palette2'][1])\n",
    "        ax2.tick_params(axis='y', labelcolor=self.color_palettes['palette2'][1])\n",
    "        ax2.set_ylim(0.7, 1.05)\n",
    "        \n",
    "        # 添加数值标签\n",
    "        for i, (week, rate) in enumerate(zip(detection_weeks, achievement_rates)):\n",
    "            ax1.text(i, week + 0.5, f'{week:.1f}周', \n",
    "                    ha='center', va='bottom', fontproperties=self.font_prop, fontsize=11)\n",
    "            ax2.text(i, rate + 0.01, f'{rate:.3f}', \n",
    "                    ha='center', va='bottom', fontproperties=self.font_prop, fontsize=11)\n",
    "        \n",
    "        plt.title('各BMI组最佳检测时点与预期达标率', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图6_检测时点优化结果.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图6：检测时点优化结果 - 已保存\")\n",
    "    \n",
    "    def plot_multifactor_correlation_heatmap(self):\n",
    "        \"\"\"图7：多因素相关性热图\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 创建相关性矩阵（基于相关性分析数据）\n",
    "        factors = self.correlation_analysis['因素'].tolist()\n",
    "        pearson_corr = self.correlation_analysis['Pearson相关系数'].tolist()\n",
    "        \n",
    "        # 构建相关性矩阵\n",
    "        n_factors = len(factors)\n",
    "        corr_matrix = np.eye(n_factors)\n",
    "        \n",
    "        # 填充相关性数据（简化处理）\n",
    "        for i in range(n_factors):\n",
    "            for j in range(n_factors):\n",
    "                if i != j:\n",
    "                    corr_matrix[i, j] = pearson_corr[i] * pearson_corr[j] * 0.5\n",
    "        \n",
    "        # 绘制热图\n",
    "        mask = np.triu(np.ones_like(corr_matrix, dtype=bool), k=1)\n",
    "        \n",
    "        sns.heatmap(corr_matrix, \n",
    "                   mask=mask,\n",
    "                   annot=True, \n",
    "                   cmap=self.create_custom_colormap('palette3'),\n",
    "                   center=0,\n",
    "                   square=True,\n",
    "                   fmt='.3f',\n",
    "                   cbar_kws={\"shrink\": 0.8})\n",
    "        \n",
    "        # 设置标签\n",
    "        factor_labels = [f.replace('BMI_最终', 'BMI').replace('多因素综合指标', '多因素指标') \n",
    "                        for f in factors]\n",
    "        plt.xticks(range(n_factors), factor_labels, \n",
    "                  fontproperties=self.font_prop, rotation=45)\n",
    "        plt.yticks(range(n_factors), factor_labels, \n",
    "                  fontproperties=self.font_prop, rotation=0)\n",
    "        \n",
    "        plt.title('多因素相关性分析热图', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图7_多因素相关性热图.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图7：多因素相关性热图 - 已保存\")\n",
    "    \n",
    "    def plot_group_characteristics_radar(self):\n",
    "        \"\"\"图8：各组特征雷达图\"\"\"\n",
    "        # 获取最佳解的数据\n",
    "        best_solution = self.detailed_results[self.detailed_results['解编号'] == 1]\n",
    "        \n",
    "        # 准备雷达图数据\n",
    "        categories = ['平均BMI', '平均年龄', '平均身高', '平均体重', '多因素指标']\n",
    "        \n",
    "        # 标准化数据到0-1范围\n",
    "        data_matrix = []\n",
    "        for _, group in best_solution.iterrows():\n",
    "            values = [\n",
    "                group['平均BMI'] / 50,  # BMI标准化\n",
    "                group['平均年龄'] / 50,  # 年龄标准化\n",
    "                group['平均身高'] / 180,  # 身高标准化\n",
    "                group['平均体重'] / 120,  # 体重标准化\n",
    "                group['多因素指标']  # 已经是0-1范围\n",
    "            ]\n",
    "            data_matrix.append(values)\n",
    "        \n",
    "        # 计算角度\n",
    "        angles = np.linspace(0, 2 * np.pi, len(categories), endpoint=False).tolist()\n",
    "        angles += angles[:1]  # 闭合\n",
    "        \n",
    "        # 为每个组创建单独的雷达图\n",
    "        colors = [self.color_palettes['palette1'][2], \n",
    "                 self.color_palettes['palette2'][2], \n",
    "                 self.color_palettes['palette3'][2]]\n",
    "        \n",
    "        for i, ((_, group), values) in enumerate(zip(best_solution.iterrows(), data_matrix)):\n",
    "            fig, ax = plt.subplots(figsize=(8, 8), subplot_kw=dict(projection='polar'))\n",
    "            \n",
    "            # 闭合数据\n",
    "            values += values[:1]\n",
    "            \n",
    "            # 绘制雷达图\n",
    "            ax.plot(angles, values, 'o-', linewidth=2.5, \n",
    "                   color=colors[i], markersize=8,\n",
    "                   markeredgecolor='white', markeredgewidth=1.5)\n",
    "            ax.fill(angles, values, alpha=0.25, color=colors[i])\n",
    "            \n",
    "            # 设置标签\n",
    "            ax.set_xticks(angles[:-1])\n",
    "            ax.set_xticklabels(categories, fontproperties=self.font_prop, fontsize=12)\n",
    "            ax.set_ylim(0, 1)\n",
    "            ax.set_yticks([0.2, 0.4, 0.6, 0.8, 1.0])\n",
    "            ax.set_yticklabels(['0.2', '0.4', '0.6', '0.8', '1.0'], fontsize=11)\n",
    "            ax.grid(True, alpha=0.3)\n",
    "            \n",
    "            plt.title(f'{group[\"组别\"]}特征分析雷达图', \n",
    "                     fontproperties=self.font_prop, fontsize=13, pad=20)\n",
    "            \n",
    "            plt.tight_layout()\n",
    "            plt.savefig(f'{self.results_dir}/图8{chr(97+i)}_{group[\"组别\"]}特征雷达图.png', \n",
    "                       dpi=300, bbox_inches='tight')\n",
    "            plt.close()\n",
    "        \n",
    "        print(\"图8a-8c：各组特征雷达图 - 已保存\")\n",
    "    \n",
    "    def plot_objective_trade_off_analysis(self):\n",
    "        \"\"\"图9：目标函数权衡分析\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 创建三个子目标的权衡分析\n",
    "        risk = self.pareto_solutions['综合风险']\n",
    "        variance = self.pareto_solutions['时间方差']\n",
    "        imbalance = self.pareto_solutions['组间不平衡']\n",
    "        \n",
    "        # 计算权衡指标\n",
    "        risk_variance_ratio = risk / (variance / 1000)  # 标准化方差\n",
    "        risk_imbalance_ratio = risk / imbalance\n",
    "        \n",
    "        # 创建散点图\n",
    "        scatter = plt.scatter(risk_variance_ratio, risk_imbalance_ratio,\n",
    "                            c=self.pareto_solutions['分组数'],\n",
    "                            cmap=self.create_custom_colormap('palette1'),\n",
    "                            s=100, alpha=0.8,\n",
    "                            edgecolors='black', linewidths=0.8)\n",
    "        \n",
    "        # 添加颜色条\n",
    "        cb = plt.colorbar(scatter, shrink=0.8)\n",
    "        cb.set_label('分组数', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        plt.xlabel('风险/方差权衡比', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('风险/不平衡权衡比', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('多目标权衡分析', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图9_多目标权衡分析.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图9：多目标权衡分析 - 已保存\")\n",
    "    \n",
    "    def plot_solution_performance_comparison(self):\n",
    "        \"\"\"图10：解决方案性能对比\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 选择前5个最佳解进行对比\n",
    "        top_solutions = self.pareto_solutions.head(5)\n",
    "        \n",
    "        solution_ids = [f'解{i}' for i in top_solutions['解编号']]\n",
    "        \n",
    "        # 标准化三个目标函数值\n",
    "        risk_norm = top_solutions['综合风险'] / top_solutions['综合风险'].max()\n",
    "        variance_norm = top_solutions['时间方差'] / top_solutions['时间方差'].max()\n",
    "        imbalance_norm = top_solutions['组间不平衡'] / top_solutions['组间不平衡'].max()\n",
    "        \n",
    "        x = np.arange(len(solution_ids))\n",
    "        width = 0.25\n",
    "        \n",
    "        # 创建分组柱状图\n",
    "        bars1 = plt.bar(x - width, risk_norm, width, \n",
    "                       color=self.color_palettes['palette1'][3], \n",
    "                       alpha=0.8, label='综合风险', \n",
    "                       edgecolor='black', linewidth=0.8)\n",
    "        \n",
    "        bars2 = plt.bar(x, variance_norm, width, \n",
    "                       color=self.color_palettes['palette2'][3], \n",
    "                       alpha=0.8, label='时间方差', \n",
    "                       edgecolor='black', linewidth=0.8)\n",
    "        \n",
    "        bars3 = plt.bar(x + width, imbalance_norm, width, \n",
    "                       color=self.color_palettes['palette3'][3], \n",
    "                       alpha=0.8, label='组间不平衡', \n",
    "                       edgecolor='black', linewidth=0.8)\n",
    "        \n",
    "        plt.xlabel('解决方案', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('标准化目标函数值', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('帕累托最优解性能对比', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.xticks(x, solution_ids, fontproperties=self.font_prop)\n",
    "        plt.legend(prop=self.font_prop, fontsize=12)\n",
    "        \n",
    "        # 添加数值标签\n",
    "        for bars in [bars1, bars2, bars3]:\n",
    "            for bar in bars:\n",
    "                height = bar.get_height()\n",
    "                plt.text(bar.get_x() + bar.get_width()/2., height + 0.01,\n",
    "                        f'{height:.3f}', ha='center', va='bottom',\n",
    "                        fontproperties=self.font_prop, fontsize=10)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图10_解决方案性能对比.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图10：解决方案性能对比 - 已保存\")\n",
    "    \n",
    "    def plot_detection_timing_distribution(self):\n",
    "        \"\"\"图11：检测时点分布分析\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 提取所有解的检测时点数据\n",
    "        all_detection_times = []\n",
    "        for _, solution in self.pareto_solutions.iterrows():\n",
    "            times_str = solution['检测时点'].strip(\"[]'\").replace(\"'\", \"\").split(\", \")\n",
    "            times = [float(t) for t in times_str if t.replace('.', '').isdigit()]\n",
    "            all_detection_times.extend([t/7 for t in times])  # 转换为周\n",
    "        \n",
    "        # 创建直方图\n",
    "        plt.hist(all_detection_times, bins=20, \n",
    "                color=self.color_palettes['palette2'][4], \n",
    "                alpha=0.8, edgecolor='black', linewidth=0.8,\n",
    "                density=True)\n",
    "        \n",
    "        # 添加核密度估计曲线\n",
    "        from scipy.stats import gaussian_kde\n",
    "        if len(all_detection_times) > 1:\n",
    "            kde = gaussian_kde(all_detection_times)\n",
    "            x_range = np.linspace(min(all_detection_times), max(all_detection_times), 100)\n",
    "            plt.plot(x_range, kde(x_range), \n",
    "                    color=self.color_palettes['palette1'][1], \n",
    "                    linewidth=3, label='核密度估计')\n",
    "        \n",
    "        plt.xlabel('检测时点（周）', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('概率密度', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('最优检测时点分布分析', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.legend(prop=self.font_prop, fontsize=12)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图11_检测时点分布分析.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图11：检测时点分布分析 - 已保存\")\n",
    "    \n",
    "    def plot_group_size_balance_analysis(self):\n",
    "        \"\"\"图12：组间样本均衡性分析\"\"\"\n",
    "        plt.figure(figsize=(8, 6))\n",
    "        \n",
    "        # 计算各解的组间不平衡程度\n",
    "        solutions = self.pareto_solutions['解编号'][:10]  # 取前10个解\n",
    "        imbalance_scores = self.pareto_solutions['组间不平衡'][:10]\n",
    "        group_numbers = self.pareto_solutions['分组数'][:10]\n",
    "        \n",
    "        # 创建散点图\n",
    "        colors = [self.color_palettes['palette3'][i % 7] for i in range(len(solutions))]\n",
    "        \n",
    "        scatter = plt.scatter(group_numbers, imbalance_scores, \n",
    "                            c=colors, s=150, alpha=0.8,\n",
    "                            edgecolors='black', linewidths=1)\n",
    "        \n",
    "        # 添加趋势线\n",
    "        z = np.polyfit(group_numbers, imbalance_scores, 1)\n",
    "        p = np.poly1d(z)\n",
    "        plt.plot(group_numbers, p(group_numbers), \n",
    "                color=self.color_palettes['palette1'][1], \n",
    "                linewidth=2.5, linestyle='--', alpha=0.8)\n",
    "        \n",
    "        plt.xlabel('分组数', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.ylabel('组间不平衡指数', fontproperties=self.font_prop, fontsize=13)\n",
    "        plt.title('分组数与组间均衡性关系', fontproperties=self.font_prop, fontsize=13)\n",
    "        \n",
    "        # 添加解编号标签\n",
    "        for i, (sol_id, x, y) in enumerate(zip(solutions, group_numbers, imbalance_scores)):\n",
    "            plt.annotate(f'解{sol_id}', (x, y), xytext=(5, 5), \n",
    "                        textcoords='offset points', \n",
    "                        fontproperties=self.font_prop, fontsize=10)\n",
    "        \n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f'{self.results_dir}/图12_组间均衡性分析.png', \n",
    "                   dpi=300, bbox_inches='tight')\n",
    "        plt.close()\n",
    "        print(\"图12：组间均衡性分析 - 已保存\")\n",
    "    \n",
    "    def generate_all_plots(self):\n",
    "        \"\"\"生成所有图片\"\"\"\n",
    "        print(\"开始生成问题三SCI风格图片...\")\n",
    "        print(\"=\"*50)\n",
    "        \n",
    "        try:\n",
    "            # 生成12类图片\n",
    "            self.plot_multifactor_bmi_hexbin()              # 图1：蜂窝图\n",
    "            self.plot_group_comparison_histograms()         # 图2a-2c：直方图对比\n",
    "            self.plot_3d_pareto_front()                     # 图3：3D帕累托前沿\n",
    "            self.plot_evolution_convergence()               # 图4：算法收敛\n",
    "            self.plot_risk_variance_scatter()               # 图5：风险方差散点\n",
    "            self.plot_detection_time_optimization()         # 图6：时点优化\n",
    "            self.plot_multifactor_correlation_heatmap()     # 图7：相关性热图\n",
    "            self.plot_group_characteristics_radar()         # 图8a-8c：雷达图\n",
    "            self.plot_objective_trade_off_analysis()        # 图9：权衡分析\n",
    "            self.plot_solution_performance_comparison()     # 图10：性能对比\n",
    "            self.plot_detection_timing_distribution()       # 图11：时点分布\n",
    "            self.plot_group_size_balance_analysis()         # 图12：均衡性分析\n",
    "            \n",
    "            print(\"=\"*50)\n",
    "            print(\"问题三SCI风格图片生成完成！\")\n",
    "            print(f\"所有图片已保存到：{self.results_dir}\")\n",
    "            \n",
    "        except Exception as e:\n",
    "            print(f\"图片生成过程中出现错误：{e}\")\n",
    "            import traceback\n",
    "            traceback.print_exc()\n",
    "\n",
    "def main():\n",
    "    \"\"\"主函数\"\"\"\n",
    "    plotter = Problem3SCIPlotter()\n",
    "    plotter.generate_all_plots()\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
